/* ======================================================================================== */
/* module aerosol.c  - functions to facilitate aerosol model selection and application      */
/*                                                                                          */
/* Description:                                                                             */
/*                                                                                          */
/* This code replaces the original set of fortran subroutines developed by M.Wang, H.Gordon,*/
/* and others (e.g., rho_a_sub_quad, linear_abc, funct_eps, load_aer, load_ss11) as well as */
/* the original get_aerosol() developed for MSl12.                                          */
/*                                                                                          */
/* The functions herein read and interpolate the aerosol model tables, which are now stored */
/* as individual HDF files per model.  Where sensor wavelengths differ from tabulated model */
/* wavelengths, interpolation is performed. Efficiencies are gained by remembering and      */
/* reusing computational results when applicable.                                           */
/*                                                                                          */
/* Primary Function:                                                                        */
/* ----------------                                                                         */
/* aerosol() - compute aerosol reflectance using specified algorithm (main interface)       */
/*                                                                                          */
/* Secondary Functions:                                                                     */
/* -------------------                                                                      */
/* wangaer() - compute aerosol reflectance using Gordon & Wang 1994 algorithm               */
/* fixedaot() - compute aerosol reflectance for fixed aot(lambda)                           */
/* fixedaer() - compute aerosol reflectance for fixed aerosol model                         */
/* get_angstrom() - compute angstrom coefficient (interface to l2_hdf_generic)              */
/* diff_tran() - compute Rayleigh-aerosol diffuse trans for selected model pair, both paths */
/*                                                                                          */
/* Supporting Functions:                                                                    */
/* --------------------                                                                     */
/* load_aermod() - loads the entire aerosol model table for the specified model list        */
/* ss_to_ms_coef() - for one model, return coefficients of function relating single         */
/*                   scattering to multiple scattering at the input geometry.               */
/* rhoas_to_rhoa() - SS aerosol reflectance to MS aerosol reflectance                       */
/* rhoa_to_rhoas() - MS aerosol reflectance to SS aerosol reflectance                       */
/* model_select_wang() - M. Wang aerosol model selection process                 .          */
/* model_select_angst() - Select model pair based on input Angstrom coefficient             */
/* model_phase() - return phase function at model wavelengths at the input geometry.        */
/* model_epsilon() - return model epsilon at input wavelengths for input geometry.          */
/* model_transmittance() - compute path Rayleigh-aerosol diffuse trans for specified model  */
/* model_taua() - compute AOT at input wavelengths using specified model                    */
/* aeroob - out-of-band water-vapor correction                                              */
/*                                                                                          */
/*                                                                                          */
/* Written By:  B. Franz, SAIC, NASA/GSFC Ocean Biology Processing Group, Summer 2004.      */
/*                                                                                          */
/* ======================================================================================== */

#include <float.h>

#include "HY1C_out.h"
#include "hy1c_util.h"

#include "l12_proto.h"
#include "config.h"

#define MAXBAND  (NBANDS+1)          /* must be >= NBANDS */
#define MAXMODEL MAXAERMOD
#define MAXSOLZ  33
#define MAXSENZ  35
#define MAXPHI   19
#define MAXSCATT 75
#define DTNTHETA 33
#define DTNBAND  MAXBAND

static float pi    = PI;
static float radeg = RADEG;
static float p0    = STDPR;

static int have_rh = 0;
static int have_sd = 0;
static int use_rh  = 0;

typedef struct aermod_struct {

    char  name[32];
    float rh;
    int   sd;

    /* angstrom exponent (nbands+1)*/
    float  *angstrom;

    /* single-scattering albedo(nbands+1), extinction coefficient(nbands+1), phase function */
    float  *albedo;
    float  *extc;
    float  **phase;

    /* quadratic coefficients for SS to MS relationship */
    float  *acost;
    float  *bcost;
    float  *ccost;

    /* Rayleigh-aerosol diffuse transmittance coeffs */
    float  **dtran_a;
    float  **dtran_b;

    /* derived quantities */
    float  **lnphase;
    float  **d2phase;

} aermodstr;

aermodstr* alloc_aermodstr( int nbands, int nscatt, int nphi, int nsolz, int nsenz, int ntheta )
{
    aermodstr  *model;

    model = (aermodstr *) malloc(sizeof(aermodstr));
    model->angstrom = (float*) malloc(nbands*sizeof(float));
    model->albedo   = (float*) malloc(nbands*sizeof(float));
    model->extc     = (float*) malloc(nbands*sizeof(float));
    model->phase    = alloc2d_float( nscatt, nbands );
    model->lnphase  = alloc2d_float( nscatt, nbands );
    model->d2phase  = alloc2d_float( nscatt, nbands );
    model->dtran_a  = alloc2d_float( ntheta, nbands );
    model->dtran_b  = alloc2d_float( ntheta, nbands );
    model->acost    = (float*) malloc( nbands*nsolz*nphi*nsenz*sizeof(float));
    model->bcost    = (float*) malloc( nbands*nsolz*nphi*nsenz*sizeof(float));
    model->ccost    = (float*) malloc( nbands*nsolz*nphi*nsenz*sizeof(float));

    return model;
}

static void free_aermodstr( aermodstr *model )
{
    free(model->extc);
    free(model->albedo);
    free(model->angstrom);
    free(model->acost);
    free(model->bcost);
    free(model->ccost);
    free2d_float(model->phase);
    free2d_float(model->lnphase);
    free2d_float(model->d2phase);
    free2d_float(model->dtran_a);
    free2d_float(model->dtran_b);
    free(model);

}

typedef struct aermodtab_struct {

    int32_t   sensorID;

    /* table dimensions */
    int32_t   nwave;
    int32_t   nmodel;
    int32_t   nsolz;
    int32_t   nsenz;
    int32_t   nphi;
    int32_t   nscatt;
    int32_t   dtran_nwave;
    int32_t   dtran_ntheta;

    /* table spectral bands and angles */
    float  *wave;
    float  *solz;
    float  *senz;
    float  *phi;
    float  *scatt;

    /* diffuse transmittance spectral bands and angles */
    float  *dtran_wave;
    float  *dtran_theta;
    float  *dtran_airmass;

    aermodstr **model;

} aermodtabstr;

static void free_aermodtabstr( aermodtabstr *aertab )
{
    int im;

    for ( im=0; im<aertab->nmodel+1; im++ )
        if ( aertab->model[im] )
            free_aermodstr(aertab->model[im]);

    free(aertab->wave);
    free(aertab->solz);
    free(aertab->senz);
    free(aertab->phi);
    free(aertab->scatt);
    free(aertab->dtran_wave);
    free(aertab->dtran_theta);
    free(aertab->dtran_airmass);
}

typedef struct alphaT_struct {
    int32_t   modnum;
    float  angstrom;
} alphaTstr;

/* aerosol table */
static aermodtabstr *aertab = NULL;

void aerosol_free( void )
{
    if ( aertab )
        free_aermodtabstr(aertab);
}


/* global variable declarations */
static int   loaded = 0;
static int   interp = 0;
static int32_t  iwatab[NBANDS];
static int32_t  iwdtab[NBANDS];

static int32_t  iwnir_s=-1;
static int32_t  iwnir_l=-1;

static int   imm50 = -1;
static int   imc50 = -1;
static int   imc70 = -1;
static int   imt90 = -1;
static int   imt99 = -1;
static int   wang_modx = 0;

static float mu0;
static float mu;
static float airmass;

static int32_t  evalmask = 0;
static int32_t  aer_opt = 0;
static float airmass_plp;
static float airmass_sph;



/* ---------------------------------------------------------------------------------------- */
/* first_deriv() - returns first derivative (dy/dx) of 1st or last array indices using a    */
/*                 4-pt Lagrangian interpolation.  NOTE: It is assumed that 4 points exist. */
/* ---------------------------------------------------------------------------------------- */
float first_deriv(float x[],float y[], int n)
{
    float a1,a2,a3,a4,a5,a6,d1;

    if (n == 0) {

        a1 = x[0]-x[1];
        a2 = x[0]-x[2];
        a3 = x[0]-x[3];
        a4 = x[1]-x[2];
        a5 = x[1]-x[3];
        a6 = x[2]-x[3];

        d1 = y[0]*(1.0/a1+1.0/a2+1.0/a3)
            - a2*a3*y[1]/(a1*a4*a5)
            + a1*a3*y[2]/(a2*a4*a6)
            - a1*a2*y[3]/(a3*a5*a6);

    } else {

        a1 = x[n-1]-x[n-4];
        a2 = x[n-1]-x[n-3];
        a3 = x[n-1]-x[n-2];
        a4 = x[n-2]-x[n-4];
        a5 = x[n-2]-x[n-3];
        a6 = x[n-3]-x[n-4];

        d1 = -a2*a3*y[n-4]/(a6*a4*a1)
            +  a1*a3*y[n-3]/(a6*a5*a2)
            -  a1*a2*y[n-2]/(a4*a5*a3)
            +  y[n-1]*(1.0/a1+1.0/a2+1.0/a3);
    }

    return(d1);
}


/* ---------------------------------------------------------------------------------------- */
/* load_aermod() - loads the entire aerosol model table for the specified model list        */
/* ---------------------------------------------------------------------------------------- */
int load_aermod(int32_t sensorID,float wave[],int32_t nwave, char models[MAXAERMOD][32], int32_t nmodels, char* configFile)
{
    char *tmp_str;
    int32 sd_id;
    int32 sds_id;
    int32 status;
    int32 rank;
    int32 nt;
    int32 dims[H4_MAX_VAR_DIMS];
    int32 nattrs;
    int32 start[4] = {0,0,0,0};

    int32_t   mwave, msolz, msenz, mphi, mscatt, dtran_mwave, dtran_mtheta;

    float d1phase1;
    float d1phaseN;
    float rh;
    int16 sd;

    char  name   [H4_MAX_NC_NAME]  = "";
    char  sdsname[H4_MAX_NC_NAME]  = "";
    char  file   [FILENAME_MAX] = "";
    char  path   [FILENAME_MAX] = "";

    int iw, im, is, iwbase, i;

    strcpy(path, configFile);
    strcat(path,"/");
    HY1C_out("\n");

    int16 sdmodel[81]={17,18,19,20,21,22,23,24,25,16,17,18,19,20,21,22,23,24,25,16,17,18,19,20,21,22,23,24,25,16,17,18,19,20,21,22,23,24,25,16,17,18,19,20,21,22,23,24,25,16,17,18,19,20,21,22,23,24,25,16,17,18,19,20,21,22,23,24,25,16,17,18,19,20,21,22,23,24,25,16,24};


    for (im=0; im<nmodels+1; im++) {

        if (im < nmodels) {
            sprintf(file, "%s%s%s_%s%s", path,
                    "aerosol/aerosol_",
                    g_cfg.aerosol_LUT,
                    models[im],
                    ".hdf");
            HY1C_out("Loading %s aerosol models from %s\n",models[im], file);

            /* Open the file and initiate the SD interface */
            sd_id = SDstart(file, DFACC_RDONLY);
            if (sd_id == -1) {
                HY1C_out("-E- %s:  Error opening file %s.\n",
                        __FILE__,file);
                exit(1);
            }

        } else {
            strcpy(file,path); strcat(file,"aerosol/aerosol_");
            strcat(file,g_cfg.aerosol_LUT); strcat(file,"_default.hdf");

            HY1C_out("Loading default aerosol model from %s\n",file);

            /* Open the file and initiate the SD interface */
            sd_id = SDstart(file, DFACC_RDONLY);
            if (sd_id == -1) {
                HY1C_out("-E- %s:  Error opening file %s.\n",
                        __FILE__,file);
                exit(1);
            }
        }

        /* read attributes which should be constant between models */

        status = SDreadattr(sd_id,SDfindattr(sd_id,"Number of Wavelengths"),
                &mwave);
        status = SDreadattr(sd_id,SDfindattr(sd_id,
                    "Number of Solar Zenith Angles"), &msolz);
        status = SDreadattr(sd_id,SDfindattr(sd_id,
                    "Number of View Zenith Angles"), &msenz);
        status = SDreadattr(sd_id,SDfindattr(sd_id,
                    "Number of Relative Azimuth Angles"), &mphi);
        status = SDreadattr(sd_id,SDfindattr(sd_id,
                    "Number of Scattering Angles"), &mscatt);
        status = SDreadattr(sd_id,SDfindattr(sd_id,
                    "Number of Diffuse Transmittance Wavelengths"  ),&dtran_mwave);
        status = SDreadattr(sd_id,SDfindattr(sd_id,
                    "Number of Diffuse Transmittance Zenith Angles"),&dtran_mtheta);

        if (im == 0) {

            int32_t   nwave = mwave;
            int32_t   nsolz = msolz;
            int32_t   nsenz = msenz;
            int32_t   nphi = mphi;
            int32_t   nscatt = mscatt;
            int32_t   dtran_nwave = dtran_mwave;
            int32_t   dtran_ntheta = dtran_mtheta;

            HY1C_out("Number of Wavelengths                          %d\n", nwave);
            HY1C_out("Number of Solar Zenith Angles                  %d\n", nsolz);
            HY1C_out("Number of View Zenith Angles                   %d\n", nsenz);
            HY1C_out("Number of Relative Azimuth Angles              %d\n", nphi);
            HY1C_out("Number of Scattering Angles                    %d\n", nscatt);

            HY1C_out("Number of Diffuse Transmittance Wavelengths    %d\n", dtran_nwave);
            HY1C_out("Number of Diffuse Transmittance Zenith Angles  %d\n", dtran_ntheta);

            // allocate the aerosol table

            if ((aertab = (aermodtabstr *)calloc(1,sizeof(aermodtabstr))) == NULL) {
                HY1C_out("Unable to allocate space for aerosol table.\n");
                exit(1);
            }

            aertab->nmodel = nmodels;
            aertab->nwave  = nwave;
            aertab->nsolz  = nsolz;
            aertab->nsenz  = nsenz;
            aertab->nphi   = nphi;
            aertab->nscatt = nscatt;

            aertab->wave  = (float *) malloc(nwave*sizeof(float));
            aertab->solz  = (float *) malloc(nsolz*sizeof(float));
            aertab->senz  = (float *) malloc(nsenz*sizeof(float));
            aertab->phi   = (float *) malloc(nphi*sizeof(float));
            aertab->scatt = (float *) malloc(nscatt*sizeof(float));

            aertab->dtran_nwave  = dtran_nwave;
            aertab->dtran_ntheta = dtran_ntheta;

            aertab->dtran_wave    = (float *) malloc(dtran_nwave*sizeof(float));
            aertab->dtran_theta   = (float *) malloc(dtran_ntheta*sizeof(float));
            aertab->dtran_airmass = (float *) malloc(dtran_ntheta*sizeof(float));

            // allocate the model tables

            if ((aertab->model = (aermodstr **)calloc(1,(nmodels+1)*sizeof(aermodstr*))) == NULL) {
                HY1C_out("Unable to allocate space for %d aerosol models.\n",nmodels+1);
                exit(1);
            }
            for (i=0; i<nmodels+1; i++) {
                if ((aertab->model[i] = alloc_aermodstr( nwave+1, nscatt, MAXPHI, MAXSOLZ, MAXSENZ, dtran_ntheta)) == NULL) {
                    HY1C_out("Unable to allocate space for aerosol model %d.\n",im);
                    exit(1);
                }
            }

            /* read SDSes which are constant between models */

            strcpy(sdsname,"wave");
            sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
            status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
            status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->wave);
            if (status != 0) {
                HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
                exit(1);
            } else {
                status = SDendaccess(sds_id);
            }

            strcpy(sdsname,"solz");
            sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
            status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
            status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->solz);
            if (status != 0) {
                HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
                exit(1);
            } else {
                status = SDendaccess(sds_id);
            }

            strcpy(sdsname,"senz");
            sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
            status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
            status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->senz);
            if (status != 0) {
                HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
                exit(1);
            } else {
                status = SDendaccess(sds_id);
            }

            strcpy(sdsname,"phi");
            sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
            status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
            status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->phi);
            if (status != 0) {
                HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
                exit(1);
            } else {
                status = SDendaccess(sds_id);
            }

            strcpy(sdsname,"scatt");
            sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
            status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
            status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->scatt);
            if (status != 0) {
                HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
                exit(1);
            } else {
                status = SDendaccess(sds_id);
            }

            strcpy(sdsname,"dtran_wave");
            sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
            status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
            status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->dtran_wave);
            if (status != 0) {
                HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
                exit(1);
            } else {
                status = SDendaccess(sds_id);
            }

            strcpy(sdsname,"dtran_theta");
            sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
            status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
            status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->dtran_theta);
            if (status != 0) {
                HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
                exit(1);
            } else {
                status = SDendaccess(sds_id);
            }

        } else {
            /*  check that all the aerosol files at least have the same
                main dimensions  */
            if( ( aertab->nwave != mwave ) || ( aertab->nsolz != msolz ) ||
                    ( aertab->nsenz != msenz ) || ( aertab->nphi != mphi ) ||
                    ( aertab->nscatt != mscatt ) ||
                    ( aertab->dtran_nwave != dtran_mwave ) ||
                    ( aertab->dtran_ntheta != dtran_mtheta ) )
            {
                HY1C_out("-E- %s, %d:  Error, Aerosol table %s\n",
                        __FILE__, __LINE__, file );
                HY1C_out( "    has different dimensions from previous tables\n"  );
                exit(1);
            }
        }

        if (im < nmodels)
            strncpy(aertab->model[im]->name,models[im],32);
        else
            strncpy(aertab->model[im]->name,"default",32);

        //20181025 debug for COCTS aertab
        //the below is origin
        //        status = SDreadattr(sd_id,SDfindattr(sd_id,"Relative Humidity"),&rh);
        //        if (status == 0) {
        //	    aertab->model[im]->rh = rh;
        //            have_rh = 1;
        //	} else
        //	    aertab->model[im]->rh = -1.0;

        //        status = SDreadattr(sd_id,SDfindattr(sd_id,"Size Distribution"),&sd);
        //        if (status == 0) {
        //	    aertab->model[im]->sd = sd;
        //            have_sd = 1;
        //	} else
        //	    aertab->model[im]->sd = -1;

        //the below is COCTS aerosol LUTs

        char rhstr[2 + 1]="0";
        strncpy(rhstr,models[im]+1,2);
        rh=atoi(rhstr);

        if (status == 0) {
            aertab->model[im]->rh = rh;
            have_rh = 1;
        } else
            aertab->model[im]->rh = -1.0;

        status = SDreadattr(sd_id,SDfindattr(sd_id,"Size Distribution"),&sd);
        // sd=sdmodel[im];
        if (status == 0) {
            aertab->model[im]->sd = sd;
            have_sd = 1;
        } else
            aertab->model[im]->sd = -1;

        strcpy(sdsname,"albedo");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
        status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->model[im]->albedo);
        if (status != 0) {
            HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
            exit(1);
        } else {
            status = SDendaccess(sds_id);
        }

        strcpy(sdsname,"extc");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
        status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->model[im]->extc);
        if (status != 0) {
            HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
            exit(1);
        } else {
            status = SDendaccess(sds_id);
        }

        /*
           strcpy(sdsname,"angstrom");
           sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
           status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
           status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->model[im]->angstrom);
           if (status != 0) {
           HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
           exit(1);
           } else {
           status = SDendaccess(sds_id);
           }
           */

        strcpy(sdsname,"phase");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
        status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->model[im]->phase[0]);
        if (status != 0) {
            HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
            exit(1);
        } else {
            status = SDendaccess(sds_id);
        }

        strcpy(sdsname,"acost");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
        status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->model[im]->acost);
        if (status != 0) {
            HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
            exit(1);
        } else {
            status = SDendaccess(sds_id);
        }

        strcpy(sdsname,"bcost");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
        status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->model[im]->bcost);
        if (status != 0) {
            HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
            exit(1);
        } else {
            status = SDendaccess(sds_id);
        }

        strcpy(sdsname,"ccost");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
        status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->model[im]->ccost);
        if (status != 0) {
            HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
            exit(1);
        } else {
            status = SDendaccess(sds_id);
        }

        strcpy(sdsname,"dtran_a");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
        status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->model[im]->dtran_a[0]);
        if (status != 0) {
            HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
            exit(1);
        } else {
            status = SDendaccess(sds_id);
        }

        strcpy(sdsname,"dtran_b");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
        status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) aertab->model[im]->dtran_b[0]);
        if (status != 0) {
            HY1C_out("-E- %s:  Error reading SDS %s from %s.\n",__FILE__,sdsname,file);
            exit(1);
        } else {
            status = SDendaccess(sds_id);
        }

        /* terminate access to the SD interface and close the file */
        status = SDend(sd_id);


        /* compute angstrom exponent for each model wavelength relative to max wavelength */
        iwbase = windex(865,aertab->wave,aertab->nwave);
        for (iw=0; iw<aertab->nwave; iw++) {
            if (iw != iwbase)
                aertab->model[im]->angstrom[iw] = -log(aertab->model[im]->extc[iw]/aertab->model[im]->extc[iwbase])/
                    log(aertab->wave[iw]/aertab->wave[iwbase]);
        }
        aertab->model[im]->angstrom[iwbase] = aertab->model[im]->angstrom[iwbase-1];

        /* precompute log of phase function and 2nd derivative (for cubic interp) */
        for (iw=0; iw<aertab->nwave; iw++) {
            for (is=0; is<aertab->nscatt; is++) {
                aertab->model[im]->lnphase[iw][is] = log(aertab->model[im]->phase[iw][is]);
            }
            d1phase1 = first_deriv(aertab->scatt,&aertab->model[im]->lnphase[iw][0],0);
            d1phaseN = first_deriv(aertab->scatt,&aertab->model[im]->lnphase[iw][0],aertab->nscatt);
            spline(aertab->scatt,
                    &aertab->model[im]->lnphase[iw][0],
                    aertab->nscatt,
                    d1phase1,
                    d1phaseN,
                    &aertab->model[im]->d2phase[iw][0]);
        }

        /* precompute airmass for diffuse transmittance */
        for (iw=0; iw<aertab->dtran_nwave; iw++) {
            for (is=0; is<aertab->dtran_ntheta; is++) {
                aertab->dtran_airmass[is] = 1.0/cos(aertab->dtran_theta[is]/radeg);
            }
        }
    }

    aertab->nmodel   = nmodels;
    aertab->sensorID = sensorID;

    /* map sensor wavelengths to table wavelengths */
    for (iw=0; iw<nwave; iw++) {
        iwatab[iw] = windex(wave[iw],aertab->wave,aertab->nwave);
        iwdtab[iw] = windex(wave[iw],aertab->dtran_wave,aertab->dtran_nwave);
        if ( abs(wave[iw]-aertab->wave[iwatab[iw]]) > 0.5 ) {
            HY1C_out("Aerosol model coefficients will be interpolated for %5.1f nm band.\n",wave[iw]);
            interp = 1;
        }
    }

    /* in case of more aertab wavelengths then sensor wavelengths */
    if (aertab->nwave != nwave)
        interp = 1;

    loaded = 1;

    /* set-up for Wang model cross-over correction */
    for (im=0; im<nmodels; im++) {
        if (strcmp(models[im],"m50") == 0) imm50 = im;
        else if (strcmp(models[im],"c50") == 0) imc50 = im;
        else if (strcmp(models[im],"c70") == 0) imc70 = im;
        else if (strcmp(models[im],"t90") == 0) imt90 = im;
        else if (strcmp(models[im],"t99") == 0) imt99 = im;
    }
    if (imm50 >= 0 && imc50 >= 0 && imc70 >= 0 && imt90 >= 0 && imt99 >= 0) {
        wang_modx = 1;
        HY1C_out("\nM. Wang model cross-over correction enabled.\n");
    }

    return(0);
}


/* ---------------------------------------------------------------------------------------- */
/* ss_to_ms_coef() - for one model, return coefficients of function relating single         */
/*                   scattering to multiple scattering at the input geometry.               */
/*                                                                                          */
/* This is effectively a C version of M. Wangs linear_a_b_c.f.  The program optimizes for   */
/* multiple calls at the same geometry by computing for all models on the first call with a */
/* new geometry.  It returns pointers to the internal static arrays of coefficients.        */
/*                                                                                          */
/* B. Franz, 1 June 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
void ss_to_ms_coef(int modnum, float solz, float senz, float phi, float **a, float **b, float **c)
{
    static float lastsolz = -999.;
    static float lastsenz = -999.;
    static float lastphi  = -999.;

    static int   computed[MAXMODEL];

    static float a_coef[MAXMODEL][MAXBAND];
    static float b_coef[MAXMODEL][MAXBAND];
    static float c_coef[MAXMODEL][MAXBAND];

    static float p, q, r;
    static float as000, as100, as010, as110, as001, as011, as101, as111;
    static float ai000, ai100, ai010, ai110, ai001, ai011, ai101, ai111;
    static float ac000, ac100, ac010, ac110, ac001, ac011, ac101, ac111;

    static int isolz1, isolz2;
    static int isenz1, isenz2;
    static int iphi1 , iphi2;

    float aphi;
    int   im, iw, i;

    /*
       if (solz >= aertab->solz[aertab->nsolz-1] || senz >= aertab->senz[aertab->nsenz-1]) {
       HY1C_out("geometry out-of-range\n");
       return(1);
       }
       */

    if (solz != lastsolz || senz != lastsenz || phi != lastphi) {

        for (im=0; im<aertab->nmodel; im++)
            computed[im] = 0;

        /* find bracketing solar indices */
        for (i=0; i<aertab->nsolz; i++) {
            if (solz < aertab->solz[i])
                break;
        }
        isolz1 = MAX(i-1,0);
        isolz2 = MIN(i,aertab->nsolz-1);
        if (isolz2 != isolz1)
            r = (solz-aertab->solz[isolz1])/(aertab->solz[isolz2]-aertab->solz[isolz1]);
        else
            r = 0.0;

        /* find bracketing view indices */
        for (i=0; i<aertab->nsenz; i++) {
            if (senz < aertab->senz[i])
                break;
        }
        isenz1 = MAX(i-1,0);
        isenz2 = MIN(i,aertab->nsenz-1);
        if (isenz2 != isenz1)
            p = (senz-aertab->senz[isenz1])/(aertab->senz[isenz2]-aertab->senz[isenz1]);
        else
            p = 0.0;

        /* find bracketing azimuth indices */
        aphi = fabs(phi);
        for (i=0; i<aertab->nphi; i++) {
            if (aphi < aertab->phi[i])
                break;
        }
        iphi1 = MAX(i-1,0);
        iphi2 = MIN(i,aertab->nphi-1);
        if (iphi2 != iphi1)
            q = (aphi-aertab->phi[iphi1])/(aertab->phi[iphi2]-aertab->phi[iphi1]);
        else
            q = 0.0;

        /* remember last geometry */
        lastsolz = solz;
        lastsenz = senz;
        lastphi  = phi;
    }

    if (!computed[modnum]) {

        im = modnum; computed[modnum] = 1;

        if (isolz2 == 0) {

            for (iw=0; iw<aertab->nwave; iw++) {

#define INDEX(iw,isol,iphi,isen) (iw*aertab->nsolz*aertab->nphi*aertab->nsenz + isol*aertab->nphi*aertab->nsenz + iphi*aertab->nsenz + isen)

                as000 = aertab->model[im]->acost[INDEX(iw,isolz1,0,isenz1)];
                as100 = aertab->model[im]->acost[INDEX(iw,isolz1,0,isenz2)];
                as001 = aertab->model[im]->acost[INDEX(iw,isolz2,0,isenz1)];
                as101 = aertab->model[im]->acost[INDEX(iw,isolz2,0,isenz2)];

                ai000 = aertab->model[im]->bcost[INDEX(iw,isolz1,0,isenz1)];
                ai100 = aertab->model[im]->bcost[INDEX(iw,isolz1,0,isenz2)];
                ai001 = aertab->model[im]->bcost[INDEX(iw,isolz2,0,isenz1)];
                ai101 = aertab->model[im]->bcost[INDEX(iw,isolz2,0,isenz2)];

                ac000 = aertab->model[im]->ccost[INDEX(iw,isolz1,0,isenz1)];
                ac100 = aertab->model[im]->ccost[INDEX(iw,isolz1,0,isenz2)];
                ac001 = aertab->model[im]->ccost[INDEX(iw,isolz2,0,isenz1)];
                ac101 = aertab->model[im]->ccost[INDEX(iw,isolz2,0,isenz2)];

                a_coef[im][iw] = (1.-p)*(1.-r)*as000 + p*r*as101
                    + (1.-p)*r*as001 + p*(1.-r)*as100;

                b_coef[im][iw] = (1.-p)*(1.-r)*ai000 + p*r*ai101
                    + (1.-p)*r*ai001 + p*(1.-q)*(1.-r)*ai100;

                c_coef[im][iw] = (1.-p)*(1.-r)*ac000 + p*r*ac101
                    + (1.-p)*r*ac001 + p*(1.-q)*(1.-r)*ac100;
            }

        } else {

            for (iw=0; iw<aertab->nwave; iw++) {

                as000 = aertab->model[im]->acost[INDEX(iw,isolz1,iphi1,isenz1)];
                as100 = aertab->model[im]->acost[INDEX(iw,isolz1,iphi1,isenz2)];
                as010 = aertab->model[im]->acost[INDEX(iw,isolz1,iphi2,isenz1)];
                as110 = aertab->model[im]->acost[INDEX(iw,isolz1,iphi2,isenz2)];
                as001 = aertab->model[im]->acost[INDEX(iw,isolz2,iphi1,isenz1)];
                as011 = aertab->model[im]->acost[INDEX(iw,isolz2,iphi2,isenz1)];
                as101 = aertab->model[im]->acost[INDEX(iw,isolz2,iphi1,isenz2)];
                as111 = aertab->model[im]->acost[INDEX(iw,isolz2,iphi2,isenz2)];

                ai000 = aertab->model[im]->bcost[INDEX(iw,isolz1,iphi1,isenz1)];
                ai100 = aertab->model[im]->bcost[INDEX(iw,isolz1,iphi1,isenz2)];
                ai010 = aertab->model[im]->bcost[INDEX(iw,isolz1,iphi2,isenz1)];
                ai110 = aertab->model[im]->bcost[INDEX(iw,isolz1,iphi2,isenz2)];
                ai001 = aertab->model[im]->bcost[INDEX(iw,isolz2,iphi1,isenz1)];
                ai011 = aertab->model[im]->bcost[INDEX(iw,isolz2,iphi2,isenz1)];
                ai101 = aertab->model[im]->bcost[INDEX(iw,isolz2,iphi1,isenz2)];
                ai111 = aertab->model[im]->bcost[INDEX(iw,isolz2,iphi2,isenz2)];

                ac000 = aertab->model[im]->ccost[INDEX(iw,isolz1,iphi1,isenz1)];
                ac100 = aertab->model[im]->ccost[INDEX(iw,isolz1,iphi1,isenz2)];
                ac010 = aertab->model[im]->ccost[INDEX(iw,isolz1,iphi2,isenz1)];
                ac110 = aertab->model[im]->ccost[INDEX(iw,isolz1,iphi2,isenz2)];
                ac001 = aertab->model[im]->ccost[INDEX(iw,isolz2,iphi1,isenz1)];
                ac011 = aertab->model[im]->ccost[INDEX(iw,isolz2,iphi2,isenz1)];
                ac101 = aertab->model[im]->ccost[INDEX(iw,isolz2,iphi1,isenz2)];
                ac111 = aertab->model[im]->ccost[INDEX(iw,isolz2,iphi2,isenz2)];

                a_coef[im][iw] = (1.-p)*(1.-q)*(1.-r)*as000 + p*q*r*as111
                    + p*(1.-q)*r*as101 + (1.-p)*q*(1.-r)*as010
                    + p*q*(1.-r)*as110 + (1.-p)*(1.-q)*r*as001
                    + (1.-p)*q*r*as011 + p*(1.-q)*(1.-r)*as100;

                b_coef[im][iw] = (1.-p)*(1.-q)*(1.-r)*ai000 + p*q*r*ai111
                    + p*(1.-q)*r*ai101 + (1.-p)*q*(1.-r)*ai010
                    + p*q*(1.-r)*ai110 + (1.-p)*(1.-q)*r*ai001
                    + (1.-p)*q*r*ai011 + p*(1.-q)*(1.-r)*ai100;

                c_coef[im][iw] = (1.-p)*(1.-q)*(1.-r)*ac000 + p*q*r*ac111
                    + p*(1.-q)*r*ac101 + (1.-p)*q*(1.-r)*ac010
                    + p*q*(1.-r)*ac110 + (1.-p)*(1.-q)*r*ac001
                    + (1.-p)*q*r*ac011 + p*(1.-q)*(1.-r)*ac100;
            }
        }
    }

    /* return pointers to coeffs for this geometry */
    *a = &a_coef[modnum][0];
    *b = &b_coef[modnum][0];
    *c = &c_coef[modnum][0];

    return;
}


/* ---------------------------------------------------------------------------------------- */
/* fresnel_coef() - computes Fresnel reflectance coefficient for specified index of refr.   */
/* ---------------------------------------------------------------------------------------- */
float fresnel_coef(float mu, float index)
{
    float sq, r2, q1;

    sq = sqrt(pow(index,2.0)-1.0+pow(mu,2.0));
    r2 = pow((mu-sq)/(mu+sq),2.0);
    q1 = (1.0-pow(mu,2.0)-mu*sq)/(1.0-pow(mu,2.0)+mu*sq);

    return(r2*(q1*q1+1.0)/2.0);
}


/* ---------------------------------------------------------------------------------------- */
/* model_phase() - return phase function at model wavelengths at the input geometry.        */
/*                                                                                          */
/* This is effectively a C version of load_ss.f by M. Wang.  The program optimizes for      */
/* multiple calls at the same geometry by computing for all models on the first call with a */
/* new geometry.  It returns a pointer to the internal static array for the requested model.*/
/*                                                                                          */
/* B. Franz, 1 June 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
float *model_phase(int modnum, float solz, float senz, float phi)
{
    static float nw = 1.334;

    static int   computed[MAXMODEL];
    static float lastsolz = -999.;
    static float lastsenz = -999.;
    static float lastphi  = -999.;
    static float phase[MAXMODEL][MAXBAND];
    static float fres1,fres2;
    static float scatt1,scatt2;

    float phase1, phase2;
    int   im, iw;

    /* recalculate only if geometry changes */

    if (solz != lastsolz || senz != lastsenz || phi != lastphi) {

        float csolz = cos(solz/radeg);
        float csenz = cos(senz/radeg);
        float cphi  = cos(phi /radeg);
        float temp;

        for (im=0; im<aertab->nmodel; im++)
            computed[im] = 0;

        /* determine scattering angles (direct and surface reflected) */
        temp   = sqrt((1.0-csenz*csenz)*(1.0-csolz*csolz)) * cphi;
        scatt1 = acos(MAX(-csenz*csolz + temp,-1.0))*radeg;
        scatt2 = acos(MIN( csenz*csolz + temp, 1.0))*radeg;

        /* compute Fresnel coefficients */
        fres1 = fresnel_coef(csenz,nw);
        fres2 = fresnel_coef(csolz,nw);

        lastsolz = solz;
        lastsenz = senz;
        lastphi  = phi;
    }

    if (!computed[modnum]) {

        im = modnum; computed[modnum] = 1;

        /* compute phase function for this geometry, all models */
        for (iw=0; iw<aertab->nwave; iw++) {
            splint(aertab->scatt,
                    &aertab->model[im]->lnphase[iw][0],
                    &aertab->model[im]->d2phase[iw][0],
                    aertab->nscatt,scatt1,&phase1);
            splint(aertab->scatt,
                    &aertab->model[im]->lnphase[iw][0],
                    &aertab->model[im]->d2phase[iw][0],
                    aertab->nscatt,scatt2,&phase2);
            //          incident diffuse   reflected   diff  dir
            phase[im][iw] = exp(phase1) + exp(phase2)*(fres1+fres2);
        }
    }

    return(&phase[modnum][0]);
}


/* ---------------------------------------------------------------------------------------- */
/* aeroob_cf() - out-of-band water-vapor scale factor                                       */
/* ---------------------------------------------------------------------------------------- */
float aeroob_cf(int modnum,float solz, float senz, float phi)
{
    static int firstCall = 1;
    static int iw1;
    static int iw2;

    float *phase;
    float  rhoas1, rhoas2;
    float  eps;
    float  cf;

    if (firstCall) {
        iw1 = windex(765,aertab->wave,aertab->nwave);
        iw2 = windex(865,aertab->wave,aertab->nwave);
        if (iw1 == iw2) iw1--;
        firstCall = 0;
    }

    phase  = model_phase(modnum,solz,senz,phi);
    rhoas1 = aertab->model[modnum]->albedo[iw1] * phase[iw1] * aertab->model[modnum]->extc[iw1];
    rhoas2 = aertab->model[modnum]->albedo[iw2] * phase[iw2] * aertab->model[modnum]->extc[iw2];
    eps    = rhoas1/rhoas2;
    cf     = log(eps)/(aertab->wave[iw2]-aertab->wave[iw1]);

    return(cf);
}


/* ---------------------------------------------------------------------------------------- */
/* aeroob - out-of-band water-vapor correction                                              */
/* ---------------------------------------------------------------------------------------- */
float aeroob(int32_t sensorID, int32_t iw, float cf, float wv)
{
    static float *a01;
    static float *a02;
    static float *a03;
    static float *a04;
    static float *a05;
    static float *a06;
    static float *a07;
    static float *a08;
    static float *a09;
    static float *a10;
    static float *a11;
    static float *a12;

    static int   firstCall = 1;

    float  f;

    if (firstCall) {
        firstCall = 0;
        rdsensorinfo(sensorID,evalmask,"oobwv01",(void **) &a01);   /* coeff #1 per sensor wave */
        rdsensorinfo(sensorID,evalmask,"oobwv02",(void **) &a02);
        rdsensorinfo(sensorID,evalmask,"oobwv03",(void **) &a03);
        rdsensorinfo(sensorID,evalmask,"oobwv04",(void **) &a04);
        rdsensorinfo(sensorID,evalmask,"oobwv05",(void **) &a05);
        rdsensorinfo(sensorID,evalmask,"oobwv06",(void **) &a06);
        rdsensorinfo(sensorID,evalmask,"oobwv07",(void **) &a07);
        rdsensorinfo(sensorID,evalmask,"oobwv08",(void **) &a08);
        rdsensorinfo(sensorID,evalmask,"oobwv09",(void **) &a09);
        rdsensorinfo(sensorID,evalmask,"oobwv10",(void **) &a10);
        rdsensorinfo(sensorID,evalmask,"oobwv11",(void **) &a11);
        rdsensorinfo(sensorID,evalmask,"oobwv12",(void **) &a12);
        HY1C_out("\nLoading water-vapor correction coefficients.\n");
    }

    f = (a01[iw] + a02[iw]*airmass + cf*(a03[iw] + a04[iw]*airmass))
        + (a05[iw] + a06[iw]*airmass + cf*(a07[iw] + a08[iw]*airmass))*wv
        + (a09[iw] + a10[iw]*airmass + cf*(a11[iw] + a12[iw]*airmass))*wv*wv;

    return(f);
}


/* ---------------------------------------------------------------------------------------- */
/* rhoa_to_rhoas() - MS aerosol reflectance to SS aerosol reflectance                       */
/*                                                                                          */
/* B. Franz, 1 June 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
int rhoa_to_rhoas(int32_t sensorID, int modnum, float solz, float senz, float phi, float wv,
        float rhoa[], float wave[], int32_t nwave, int iw1, int iw2, float rhoas[])
{
    float *ac,*bc,*cc;
    double a,b,c;
    double f;
    int   iw;
    int   iwtab;
    float cf;
    int status = 0;

    ss_to_ms_coef(modnum,solz,senz,phi,&ac,&bc,&cc);

    cf = aeroob_cf(modnum,solz,senz,phi);

    for (iw=iw1; iw<=iw2; iw++) {
        if (rhoa[iw] < 1.e-20)
            rhoas[iw] = rhoa[iw];
        else {
            iwtab = iwatab[iw];
            a = (double) ac[iwtab];
            b = (double) bc[iwtab];
            c = (double) cc[iwtab];
            f = b*b - 4*c*( a - log((double)rhoa[iw]) );
            if (f > 0.0) {
                if (fabs(c) > 1.e-20) {
                    rhoas[iw] = exp(0.5*(-b+sqrt(f))/c);
                } else if (fabs(a) > 1.e-20 && fabs(b) > 1.e-20) {
                    rhoas[iw] = pow(rhoa[iw]/a, 1./b);
                } else {
                    status = 1;
                    break;
                }
                rhoas[iw] = rhoas[iw]/aeroob(sensorID,iw,cf,wv);
                if (!finite(rhoas[iw]) || rhoas[iw] < 1.e-20) {
                    status=1;
                    break;
                }
            } else {
                status=1;
                break;
            }
        }
    }

    // return input values and failure status if any wavelengths failed

    if (status != 0) {
        for (iw=iw1; iw<=iw2; iw++)
            rhoas[iw] = rhoa[iw];
    }

    return(status);
}



/* ---------------------------------------------------------------------------------------- */
/* rhoas_to_rhoa() - SS aerosol reflectance to MS aerosol reflectance                       */
/*                                                                                          */
/* B. Franz, 1 June 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
void rhoas_to_rhoa(int32_t sensorID, int modnum, float solz, float senz, float phi, float wv,
        float rhoas[], float wave[], int32_t nwave, int iw1, int iw2, float rhoa[])
{
    float *ac,*bc,*cc;
    float a,b,c;
    float lnrhoas;
    int   iw;
    int   iwtab;
    float cf;

    ss_to_ms_coef(modnum,solz,senz,phi,&ac,&bc,&cc);

    cf = aeroob_cf(modnum,solz,senz,phi);

    for (iw=iw1; iw<=iw2; iw++) {


        /* these changes ensure that rhoa1 -> rhoas -> rhoa2 == rhoa1 */
        /* but, that changes everything, slightly (tau)               */

        if (rhoas[iw] < 1.e-20)
            rhoa[iw] = rhoas[iw];
        else {
            iwtab = iwatab[iw];
            a = ac[iwtab];
            b = bc[iwtab];
            c = cc[iwtab];
            lnrhoas = log(rhoas[iw]*aeroob(sensorID,iw,cf,wv));
            rhoa[iw] = exp(a + b*lnrhoas + c*lnrhoas*lnrhoas);
        }

        /*
           iwtab = iwatab[iw];
           a = ac[iwtab];
           b = bc[iwtab];
           c = cc[iwtab];
           lnrhoas = log(rhoas[iw]);
           rhoa[iw] = exp(a + b*lnrhoas + c*lnrhoas*lnrhoas)
         * aeroob(sensorID,iw,cf,wv);
         */
    }

    return;
}


/* ---------------------------------------------------------------------------------------- */
/* model_epsilon() - return model epsilon at input wavelengths for input geometry.          */
/*                                                                                          */
/* If the input wavelengths are not equal to the model wavelengths, the model epsilon will  */
/* be interpolated to the input wavelengths.  It is assumed that the longest (last) input   */
/* wavelength is equivalent to the longest wavelength of the model table.  Hence,           */
/* the function should always be called with the full array of input sensor wavelengths.    */
/*                                                                                          */
/* The program optimizes for multiple calls at the same geometry by computing for all       */
/* models on the first call with a new geometry.  It returns a pointer to the internal      */
/* static arrays of epsilon for the requested model.                               .        */
/*                                                                                          */
/* B. Franz, 1 June 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
float *model_epsilon(int modnum, int32_t iwnir_l, float wave[], int32_t nwave, float solz, float senz, float phi)
{
    static float lastsolz = -999.;
    static float lastsenz = -999.;
    static float lastphi  = -999.;
    static int32_t  lastiwl  = -999;
    static int   lastmod  = -999;
    static float epsilon[MAXMODEL][MAXBAND];

    /* recalculate only if geometry changes */

    if (modnum != lastmod || solz != lastsolz || senz != lastsenz || phi != lastphi || iwnir_l != lastiwl) {

        int    iwnir  = iwatab[iwnir_l];
        float *phase;
        float  lneps[MAXBAND];
        float  rhoas1, rhoas2;
        int    im, iw, iwtab;

        im = modnum;
        phase = model_phase(im,solz,senz,phi);
        for (iw=0; iw<aertab->nwave; iw++) {
            rhoas1 = aertab->model[im]->albedo[iw]    * phase[iw]    * aertab->model[im]->extc[iw];
            rhoas2 = aertab->model[im]->albedo[iwnir] * phase[iwnir] * aertab->model[im]->extc[iwnir];
            epsilon[im][iw] = rhoas1/rhoas2;
            if (interp)
                lneps[iw] = log(epsilon[im][iw]);
        }
        if (interp) {
            for (iw=0; iw<nwave; iw++) {
                iwtab = iwatab[iw];
                if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0)
                    epsilon[im][iw] = exp(linterp(aertab->wave,lneps,aertab->nwave,wave[iw]));
                else
                    epsilon[im][iw] = exp(lneps[iwtab]);
            }
        }

        lastsolz = solz;
        lastsenz = senz;
        lastphi  = phi;
        lastiwl  = iwnir_l;
        lastmod  = modnum;
    }
    return(&epsilon[modnum][0]);
}


/* ---------------------------------------------------------------------------------------- */
/* model_select_wang() - M. Wang aerosol model selection process                 .          */
/*                                                                                          */
/* B. Franz, 1 June 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
int model_select_wang(int32_t sensorID, float wave[], int32_t nwave, float solz, float senz, float phi,
        float wv, float rhoa[], int32_t iwnir_s, int32_t iwnir_l,
        int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir)
{
    float rhoas [NBANDS];
    float eps_ret  [MAXMODEL];
    float eps_mod  [MAXMODEL];
    float eps_err  [MAXMODEL];
    int   imod     [MAXMODEL];
    int   nmod = aertab->nmodel;
    float eps_ave;
    float *eps;
    float err_m;
    float err_p;
    int   im, iim;
    int   eps_flg = 0;
    float wt;
    float tot_err;
    int   itmp;


    /* get theoretical and retrieved epsilon for each model, and */
    /* compute average retrieved epsilon                         */

    eps_ave = 0.0;
    for (im=0; im<nmod; im++) {

        /* retrieve epsilon in NIR assuming each model */
        rhoa_to_rhoas(sensorID,im,solz,senz,phi,wv,rhoa,wave,nwave,iwnir_s,iwnir_l,rhoas);
        eps_ret[im] = rhoas[iwnir_s]/rhoas[iwnir_l];

        /* get model epsilon for each model at this geometry */
        eps = model_epsilon(im,iwnir_l,wave,nwave,solz,senz,phi);
        eps_mod[im] = eps[iwnir_s];

        eps_ave += eps_ret[im];
    }
    eps_ave /= nmod;

    /* determine average retrieved epsilon for the four retrievals which most */
    /* closely match the theoretical model epsilons. the model set is reduced */
    /* from the full suite to the final four by iterative outlier rejection.  */

    while (nmod > 4) {

        /* compute differences between retrieved and model epsilon */
        for (im=0; im<aertab->nmodel; im++) {
            imod[im]    = im;
            eps_err[im] = eps_ave - eps_mod[im];
        }

        /* sort model indices by smallest to largest absolute differences */
        for (im=0; im<aertab->nmodel-1; im++)
            for (iim=im+1; iim<aertab->nmodel; iim++)
                if ( fabs(eps_err[imod[im]]) > fabs(eps_err[imod[iim]]) ) {
                    itmp      = imod[im ];
                    imod[im ] = imod[iim];
                    imod[iim] = itmp;
                }

        /* recompute average retrieved epsilon over the n-2 closest models  */
        /* averaging is done as a weighted mean with wt=1-|eps_err|/tot_err */
        /* note that the sum of the weights is equal to n-1                 */

        nmod = nmod - 2;

        tot_err = 0.0;
        for (iim=0; iim<nmod; iim++) {
            im = imod[iim];
            tot_err += fabs(eps_err[im]);
        }

        eps_ave = 0.0;
        for (iim=0; iim<nmod; iim++) {
            im = imod[iim];
            wt = 1.0-fabs(eps_err[im])/tot_err;
            eps_ave += eps_ret[im]*wt;
        }
        eps_ave /= (nmod-1);
    }

    /* now select the two models which bracket eps_ave  */
    err_m   = -999.;
    err_p   =  999.;
    *modmin = -1;
    *modmax = -1;
    for (im=0; im<aertab->nmodel; im++) {
        eps_err[im] = eps_ave - eps_mod[im];
        if (eps_err[im] >= 0.0) {
            if (eps_err[im] < err_p) {
                err_p   = eps_err[im];
                *modmin = im;
            }
        } else {
            if (eps_err[im] > err_m) {
                err_m   = eps_err[im];
                *modmax = im;
            }
        }
    }

    /* M. Wang's model cross-over correction */
    if (wang_modx && eps_mod[imc50] > eps_mod[imt99]) {
        if (*modmax == imt90)
            *modmin = imt99;
        else if (*modmax == imc50 && *modmin == imt99)
            *modmin = imc70;
        else if (*modmin == imm50)
            *modmax = imc50;
    }


    /* compute interpolation ratio */
    if (*modmin < 0) {
        /* no lower-bounding model */
        *modmin = *modmax;
        *modrat = 0.0;
        eps_flg = -1;
    } else if (*modmax < 0) {
        /* no upper-bounding model */
        *modmax = *modmin;
        *modrat = 0.0;
        eps_flg =  1;
    } else
        *modrat = (eps_ave-eps_mod[*modmin])/(eps_mod[*modmax]-eps_mod[*modmin]);

    /* return retrieved epsilon */
    *epsnir = eps_ave;

    return(eps_flg);
}

/* ---------------------------------------------------------------------------------------- */
/* model_select_wang2() - M. Wang aerosol model selection process                 .          */
/*                                                                                          */
/* B. Franz, 1 June 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
int model_select_wang2(int32_t sensorID, float wave[], int32_t nwave, int32_t nmodel, int32_t mindx[],
        float solz, float senz, float phi,
        float wv, float rhoa[], int32_t iwnir_s, int32_t iwnir_l,
        int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir)
{
    float rhoas [NBANDS];
    float eps_ret  [MAXMODEL];
    float eps_mod  [MAXMODEL];
    float eps_err  [MAXMODEL];
    int   imod     [MAXMODEL];
    int   nmod = nmodel;
    float eps_ave;
    float *eps;
    float err_m;
    float err_p;
    int   jm, im, iim;
    int   eps_flg = 0;
    float wt;
    float tot_err;
    int   itmp;

    *modmin =  -1;
    *modmax =  -1;
    *modrat = 0.0;
    *epsnir = 0.0;

    /* get theoretical and retrieved epsilon for each model, and */
    /* compute average retrieved epsilon                         */

    eps_ave = 0.0;
    for (jm=0; jm<nmod; jm++) {

        im = mindx[jm];

        /* retrieve epsilon in NIR assuming each model */
        rhoa_to_rhoas(sensorID,im,solz,senz,phi,wv,rhoa,wave,nwave,iwnir_s,iwnir_l,rhoas);
        eps_ret[jm] = rhoas[iwnir_s]/rhoas[iwnir_l];

        /* get model epsilon for each model at this geometry */
        eps = model_epsilon(im,iwnir_l,wave,nwave,solz,senz,phi);
        eps_mod[jm] = eps[iwnir_s];

        eps_ave += eps_ret[jm];
    }
    eps_ave /= nmod;


    /* determine average retrieved epsilon for the four retrievals which most */
    /* closely match the theoretical model epsilons. the model set is reduced */
    /* from the full suite to the final four by iterative outlier rejection.  */

    while (nmod > 4) {

        /* compute differences between retrieved and model epsilon */
        for (im=0; im<nmodel; im++) {
            imod[im]    = im;
            eps_err[im] = eps_ave - eps_mod[im];
        }

        /* sort model indices by smallest to largest absolute differences */
        for (im=0; im<nmodel-1; im++) {
            for (iim=im+1; iim<nmodel; iim++)
                if ( fabs(eps_err[imod[im]]) > fabs(eps_err[imod[iim]]) ) {
                    itmp      = imod[im ];
                    imod[im ] = imod[iim];
                    imod[iim] = itmp;
                }
        }

        /* recompute average retrieved epsilon over the n-2 closest models  */
        /* averaging is done as a weighted mean with wt=1-|eps_err|/tot_err */
        /* note that the sum of the weights is equal to n-1                 */

        nmod = nmod - 2;

        tot_err = 0.0;
        for (iim=0; iim<nmod; iim++) {
            im = imod[iim];
            tot_err += fabs(eps_err[im]);
        }

        eps_ave = 0.0;
        for (iim=0; iim<nmod; iim++) {
            im = imod[iim];
            wt = 1.0-fabs(eps_err[im])/tot_err;
            eps_ave += eps_ret[im]*wt;
        }
        eps_ave /= (nmod-1);
    }

    /* now select the two models which bracket eps_ave  */
    err_m   = 0-FLT_MAX;
    err_p   = FLT_MAX;
    for (im=0; im<nmodel; im++) {
        eps_err[im] = eps_ave - eps_mod[im];
        if (eps_err[im] >= 0.0) {
            if (eps_err[im] < err_p) {
                err_p   = eps_err[im];
                *modmin = im;
            }
        } else {
            if (eps_err[im] > err_m) {
                err_m   = eps_err[im];
                *modmax = im;
            }
        }
    }

    /* M. Wang's model cross-over correction */
    if (wang_modx && eps_mod[imc50] > eps_mod[imt99]) {
        if (*modmax == imt90)
            *modmin = imt99;
        else if (*modmax == imc50 && *modmin == imt99)
            *modmin = imc70;
        else if (*modmin == imm50)
            *modmax = imc50;
    }

    /* compute interpolation ratio */
    if (*modmin < 0) {
        /* no lower-bounding model */
        *modmin = *modmax;
        *modrat = 0.0;
        eps_flg = -1;
    } else if (*modmax < 0) {
        /* no upper-bounding model */
        *modmax = *modmin;
        *modrat = 0.0;
        eps_flg =  1;
    } else
        *modrat = (eps_ave-eps_mod[*modmin])/(eps_mod[*modmax]-eps_mod[*modmin]);

    *modmin = mindx[*modmin];
    *modmax = mindx[*modmax];

    /* return retrieved epsilon */
    *epsnir = eps_ave;

    return(eps_flg);
}


/* ---------------------------------------------------------------------------------------- */
/* model_select_angst() - Select model pair based on input Angstrom coefficient             */
/*                                                                                          */
/* B. Franz, 1 August 2004.                                                                 */
/* ---------------------------------------------------------------------------------------- */
int compalphaT(alphaTstr *x, alphaTstr *y) {return(x->angstrom < y->angstrom ? -1 : 1);}

void model_select_angstrom(float angstrom, int32_t *modmin, int32_t *modmax, float *modrat)
{
    static alphaTstr alphaT[MAXAERMOD];
    static int firstCall = 1;

    int im, im1, im2;

    if (firstCall) {

        int ib = windex(520,aertab->wave,aertab->nwave);

        /* grab angstrom coefficients and sort in ascending order */
        for (im=0; im<aertab->nmodel; im++) {
            alphaT[im].angstrom = aertab->model[im]->angstrom[ib];
            alphaT[im].modnum = im;
        }
        qsort(alphaT,aertab->nmodel,sizeof(alphaTstr),
                (int (*)(const void *,const void *)) compalphaT);

        firstCall = 0;
    }

    for (im=0; im<aertab->nmodel; im++) {
        if (angstrom < alphaT[im].angstrom)
            break;
    }
    im1 = MAX(MIN(im-1,aertab->nmodel-1),0);
    im2 = MAX(MIN(im  ,aertab->nmodel-1),0);

    *modmin = alphaT[im1].modnum;
    *modmax = alphaT[im2].modnum;

    if (im1 == im2)
        *modrat = 1.0;
    else
        *modrat = (angstrom-alphaT[im1].angstrom) /
            (alphaT[im2].angstrom-alphaT[im1].angstrom);

    return;

}

static int order_models(const void *p1, const void *p2);
static void print_rhaer(void)
{
    static float rhtab[MAXAERMOD];
    static int   sdtab[MAXAERMOD];
    float lastrh = -1.0;
    int   lastsd = -1;
    int   irh, isd;

    if (!have_rh || !have_sd) {
        HY1C_out("-E- %s line %d: This aerosol selection method requires models with a Relative Humidity attribute.\n",
                __FILE__,__LINE__);
        //exit(1);
    }
    // need in order of rh and sd within rh
    qsort(aertab->model,aertab->nmodel,sizeof(aermodstr*),(int (*)(const void *,const void *)) order_models);
    int nsd = 0;
    int nrh = 0;
    int im = 0;
    for (im=0; im<aertab->nmodel; im++) {
        if (aertab->model[im]->rh != lastrh) {
            rhtab[nrh] = aertab->model[im]->rh;
            lastrh     = rhtab[nrh];
            nrh++;
        }
        if (nrh == 1 && aertab->model[im]->sd != lastsd) {
            sdtab[nsd] = aertab->model[im]->sd;
            lastsd     = sdtab[nsd];
            nsd++;
        }
    }
    if (nrh*nsd != aertab->nmodel) {
        HY1C_out("-E- %s line %d: number of humidities (%d) x number of size distributions (%d) must equal number of models (%d).\n",
                __FILE__,__LINE__,nrh,nsd,aertab->nmodel);
        //exit(1);
    } else {
        HY1C_out("%d aerosol models: %d humidities x %d size fractions\n",aertab->nmodel,nrh,nsd);
        for (irh=0; irh<nrh; irh++) {
            for (isd=0; isd<nsd; isd++) {
                im = irh*nsd+isd;
                HY1C_out("model %d, rh=%f, sd=%d, alpha=%f, name=%s\n",
                        im,aertab->model[im]->rh, aertab->model[im]->sd, aertab->model[im]->angstrom[1],aertab->model[im]->name);
            }
        }
    }
}

int HY1C_cal_La(l2str *l2rec_1, l2str *l2rec_2,
        float wave[], int32_t nwave,
        char models[MAXAERMOD][32],
        int32_t nmodels, char* configFile)
{
    l2str *l2rec_H = l2rec_1;
    l2str *l2rec_M = l2rec_2;

    /* calculate aerindex vvvvvv */
    int modnum = 0;

    if ((g_cfg.DPG_method == DPG_TRBMC) &&
            (l2rec_H->field_aer_modnum != NULL) &&
            (*(l2rec_H->field_aer_modnum) != -999)) {
        modnum = *(l2rec_H->field_aer_modnum) - g_cfg.model_number_base;
    } else {
        modnum = *(l2rec_M->aermodmin) - g_cfg.model_number_base;
    }

    if (modnum < 0 || modnum == BAD_INT)
        modnum = nmodels - 1; // use the last one and the default one.

    *(l2rec_H->aerindex) = modnum;
    /* calculate aerindex ^^^^^^ */

    float solz = *(l2rec_H->solz);
    float senz = *(l2rec_H->senz);
    float wv = *(l2rec_H->wv);
    float *Fo = l2rec_H->Fo;
    int32_t sensorID = l2rec_H->sensorID;

    /* load aerosol model tables */
    if (!loaded) {
        int32_t im;
        load_aermod(sensorID, wave, nwave, models,
                nmodels, configFile);
        //for (im=0; im<aertab->nmodel; im++) mindx[im] = im;
        if (have_rh && aertab->nmodel >= 30) {
            HY1C_out("Limiting aerosol models based on RH.\n");
            use_rh = 1;
            print_rhaer();
        }
    }

    if (g_cfg.DPG_method == DPG_TRBMC)
        l2rec_H->aot_new = l2rec_H->field_AOT;
    else if (g_cfg.DPG_method == DPG_TRBMCC) {

        float *aot = l2rec_M->taua;

        switch (g_cfg.sensor_id) {
        case SENSOR_COCTS:
            line_interpolation_float(g_cfg.wave_hy1c, g_cfg.hy1c_nwave,
                    g_cfg.wave_modis, g_cfg.modis_nwave,
                    aot, l2rec_H->aot_new);
            break;
        case SENSOR_CZI:
        case SENSOR_UVI: {
            int hy1c_nwave = g_cfg.hy1c_nwave;
            int modis_nwave = g_cfg.modis_nwave;

            int iwnir = iwatab[hy1c_nwave - 1];
            if (g_cfg.sensor == SENSOR_UVI)
                iwnir = iwatab[2 - 1];

            for (int iw = 0; iw < nwave; iw++) {
                l2rec_H->aot_new[iw] = aot[modis_nwave - 1] *
                    aertab->model[modnum]->extc[iwatab[iw]] /
                    aertab->model[modnum]->extc[iwnir];
            }
            break;
        }
        }
    }

    float rhoas[8] = {0};
    float rhoa[8] = {0};
    float phi = *(l2rec_H->delphi);
    float *phase  = model_phase(modnum, solz, senz, phi);
    float  csolz  = cos(solz/radeg);
    float  csenz  = cos(senz/radeg);

    l2rec_H->phase = phase;

    for (int iw = 0; iw < nwave; iw++) {
        rhoas[iw] = l2rec_H->aot_new[iw] * phase[iw] * aertab->model[modnum]->albedo[iwatab[iw]] /(4.0 * csolz * csenz);
    }

    rhoas_to_rhoa(sensorID, modnum, solz, senz, phi, wv, rhoas, wave, nwave, 0, nwave - 1, rhoa);

    for (int iw = 0; iw < nwave; iw++) {
        (l2rec_H->La)[iw] = rhoa[iw] * Fo[iw] * cos(solz / radeg) / PI;
    }

    return 0;
}

int HY1C_cal_La2(float gain[], size_t length,
                 l2str *l2rec_1, l2str *l2rec_2,
                 aestr *aerec)
{
    float gain_s = BAD_FLT;
    float gain_l = BAD_FLT;

    if (g_cfg.again.gain[0] == BAD_FLT && g_cfg.again.gain[1] == BAD_FLT) {
        float gain_s = gain[length - 2]; /* 750 */
        float gain_l = gain[length - 1]; /* 865 */
        LOG_DBG("use generated gain %f %f\n", gain_s, gain_l);
    } else {
        float gain_s = g_cfg.again.gain[0]; /* 750 */
        float gain_l = g_cfg.again.gain[1]; /* 865 */
        LOG_DBG("use gain from input %f %f\n", gain_s, gain_l);
    }

    l2rec_1->Lt[length - 2] *= gain_s;
    l2rec_1->Lt[length - 1] *= gain_l;

    return atmocor2(l2rec_1, aerec, 0);
}

int HY1C_cal_gain(l2str *l2rec_1, l2str *l2rec_2,
        float wave[], int32_t nwave)
{
    l2str *l2rec_H = l2rec_1;
    l2str *l2rec_M = l2rec_2;

    float *Lt_theroy = l2rec_H->Lt_theroy;

    int iw = 0;

    if (g_cfg.DPG_method == DPG_TRBMC) {
        l2rec_H->Lw = l2rec_H->field_Lw;
    } else {
        float *Lw = l2rec_M->Lw;
        line_interpolation_float(wave, nwave, g_cfg.wave_modis, g_cfg.modis_nwave, Lw, l2rec_H->Lw);
        for (int i = 0; i < nwave; i++) {
            if (l2rec_H->Lw[i] < 0)
                l2rec_H->Lw[i] = 0;
        }
    }

    for (iw = 0; iw < nwave; iw++) {
        float tLw = l2rec_H->Lw[iw] * (l2rec_H->t_sol[iw]) * (l2rec_H->t_sen[iw]) *
            (l2rec_H->tg_sol[iw]) * (l2rec_H->tg_sen[iw]);

        Lt_theroy[iw] = ((l2rec_H->Lr)[iw] + (l2rec_H->La)[iw] +
                (l2rec_H->TLg)[iw] + (l2rec_H->tLf)[iw]) *
            (l2rec_H->tg_sol[iw]) * (l2rec_H->tg_sen[iw])
            + tLw;

        l2rec_H->gain[iw] = Lt_theroy[iw] / l2rec_H->Lt[iw];
    }
    return 0;
}

int HY1C_cal_gain2(l2str *l2rec_1, l2str *l2rec_2,
        float wave[], int32_t nwave)
{
    l2str *l2rec_H = l2rec_1;
    l2str *l2rec_M = l2rec_2;

    float *Lt_theroy = l2rec_H->Lt_theroy;

    int iw = 0;

    if (g_cfg.DPG_method == DPG_TRBMC) {
        l2rec_H->Lw = l2rec_H->field_Lw;
    } else {
        float *Lw = l2rec_M->Lw;
        line_interpolation_float(wave, nwave, g_cfg.wave_modis, g_cfg.modis_nwave, Lw, l2rec_H->Lw);
        for (int i = 0; i < nwave; i++) {
            if (l2rec_H->Lw[i] < 0)
                l2rec_H->Lw[i] = 0;
        }
    }

    convl21(l2rec_H, NULL, 0, 0, &(g_cfg.input), Lt_theroy, NULL);

    for (iw = 0; iw < nwave; iw++) {
        if (Lt_theroy[iw] == BAD_FLT) {
            l2rec_H->gain[iw] = BAD_FLT;
            continue;
        }

        l2rec_H->gain[iw] = Lt_theroy[iw] / l2rec_H->Lt[iw];
    }
    return 0;
}

/* ---------------------------------------------------------------------------------------- */
/* model_taua() - compute AOT at input wavelengths using specified model                    */
/* ---------------------------------------------------------------------------------------- */
void model_taua(int32_t sensorID, int modnum, float wave[], int32_t nwave, int32_t iwnir_l, float rhoa[],
        float solz, float senz, float phi, float wv, float taua[])
{
    float  aot  [MAXBAND];
    float  lnaot[MAXBAND];
    float  rhoas[NBANDS];

    int    iwnir  = iwatab[iwnir_l];
    float *phase  = model_phase(modnum,solz,senz,phi);
    float  csolz  = cos(solz/radeg);
    float  csenz  = cos(senz/radeg);
    int    iw, iwtab;

    /* get SS aerosol reflectance at longest sensor wavelength */
    rhoa_to_rhoas(sensorID,modnum,solz,senz,phi,wv,rhoa,wave,nwave,iwnir_l,iwnir_l,rhoas);

    /* get aerosol optical thickness at longest sensor wavelength */
    aot[iwnir] = rhoas[iwnir_l]*(4.0*csolz*csenz)/(phase[iwnir]*aertab->model[modnum]->albedo[iwnir]);

    /* get aerosol optical thickness at all other table wavelengths */
    for (iw=0; iw<aertab->nwave; iw++) {
        /* note to self: i actually have aot(865) for czcs at this point */
        aot[iw] = aot[iwnir] * aertab->model[modnum]->extc[iw] / aertab->model[modnum]->extc[iwnir];
        if (interp)
            lnaot[iw] = log(aot[iw]);
    }

    /* interpolate aot from table to sensor wavelengths */
    if (interp) {
        for (iw=0; iw<nwave; iw++) {
            iwtab = iwatab[iw];
            if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0)
                taua[iw] = exp(linterp(aertab->wave,lnaot,aertab->nwave,wave[iw]));
            else
                taua[iw] = aot[iwtab];
        }
    } else
        for (iw=0; iw<nwave; iw++)
            taua[iw] = aot[iw];

    return;
}


/* ---------------------------------------------------------------------------------------- */
/* model_transmittance() - compute path Rayleigh-aerosol diffuse trans for specified model  */
/* ---------------------------------------------------------------------------------------- */
void model_transmittance(int modnum, float wave[], int32_t nwave,
        float theta, float taua[], float dtran[])
{
    static int firstCall = 1;
    static float *intexp;
    static int   *inttst;

    int   i1, i2, i;
    int   iw, iwtab;
    float a1, b1;
    float a2, b2;
    float x1,x2;
    float y1,y2;
    float xbar;
    float wt;

    if (firstCall) {
        float taur1, um1;
        float taur2, um2;
        firstCall = 0;

        intexp = (float *) malloc(nwave*sizeof(float));
        inttst =   (int *) malloc(nwave*sizeof(int));

        for (iw=0; iw<nwave; iw++) {
            intexp[iw] = 1.0;
            inttst[iw] = 0;
            iwtab = iwdtab[iw];
            if (fabs(wave[iw] - aertab->dtran_wave[iwtab]) > 0.51) {
                um1 = aertab->dtran_wave[iwtab]/1000.0;
                um2 = wave[iw]/1000.0;
                taur1 = 0.008569*pow(um1,-4)*(1.0+(0.0113*pow(um1,-2))+(0.00013*pow(um1,-4)));
                taur2 = 0.008569*pow(um2,-4)*(1.0+(0.0113*pow(um2,-2))+(0.00013*pow(um2,-4)));
                intexp[iw] = taur2/taur1;
                inttst[iw] = 1;
                HY1C_out("Interpolating diffuse transmittance for %d from %f by %f\n",
                        (int) wave[iw],aertab->dtran_wave[iwtab],intexp[iw]);
            }
        }
    }

    /* find bracketing zenith angle indices */
    for (i=0; i<aertab->dtran_ntheta; i++) {
        if (theta < aertab->dtran_theta[i])
            break;
    }
    if (i == aertab->dtran_ntheta) {
        i1 = i-1;
        i2 = i-1;
        wt = 0.0;
    } else {
        i1 = MIN(MAX(i-1,0),aertab->dtran_ntheta-2);
        i2 = i1+1;
        x1   = aertab->dtran_airmass[i1];
        x2   = aertab->dtran_airmass[i2];
        xbar = 1.0/cos(theta/radeg);
        wt   = (xbar-x1)/(x2-x1);
    }

    /* use coefficients of nearest wavelength */
    for (iw=0; iw<nwave; iw++) {

        iwtab = iwdtab[iw];

        a1  = aertab->model[modnum]->dtran_a[iwtab][i1];
        b1  = aertab->model[modnum]->dtran_b[iwtab][i1];

        a2  = aertab->model[modnum]->dtran_a[iwtab][i2];
        b2  = aertab->model[modnum]->dtran_b[iwtab][i2];

        if (inttst[iw]) {
            a1 = pow(a1,intexp[iw]);
            a2 = pow(a2,intexp[iw]);
        }

        y1 = a1 * exp(-b1*taua[iw]);
        y2 = a2 * exp(-b2*taua[iw]);

        dtran[iw] = MAX(MIN( (1.0-wt)*y1 + wt*y2, 1.0), 1e-5);
    }

    return;
}


/* ---------------------------------------------------------------------------------------- */
/* diff_tran() - compute Rayleigh-aerosol diffuse trans for selected model pair, both paths */
/* ---------------------------------------------------------------------------------------- */
void diff_tran(int32_t sensorID,float wave[], int32_t nwave, int32_t iwnir_l, float solz, float senz, float phi,
        float wv, float pr, float taur[], int32_t modmin, int32_t modmax, float modrat, float rhoa[],
        float taua[], float tsol[], float tsen[])
{
    int iw;
    float tauamin[NBANDS];
    float tauamax[NBANDS];
    float tsolmin[NBANDS];
    float tsolmax[NBANDS];
    float tsenmin[NBANDS];
    float tsenmax[NBANDS];

    /* get AOT per band for each model */
    model_taua(sensorID,modmin,wave,nwave,iwnir_l,rhoa,solz,senz,phi,wv,tauamin);
    model_taua(sensorID,modmax,wave,nwave,iwnir_l,rhoa,solz,senz,phi,wv,tauamax);

    /* get diff trans sun to ground, per band for each model */
    model_transmittance(modmin,wave,nwave,solz,tauamin,tsolmin);
    model_transmittance(modmax,wave,nwave,solz,tauamax,tsolmax);

    /* get diff trans ground to sensor, per band for each model */
    model_transmittance(modmin,wave,nwave,senz,tauamin,tsenmin);
    model_transmittance(modmax,wave,nwave,senz,tauamax,tsenmax);

    /* interpolate and pressure correct */
    for (iw=0; iw<nwave; iw++) {
        taua[iw] = tauamin[iw]*(1.0-modrat) + tauamax[iw]*modrat;
        tsol[iw] = tsolmin[iw]*(1.0-modrat) + tsolmax[iw]*modrat;
        tsen[iw] = tsenmin[iw]*(1.0-modrat) + tsenmax[iw]*modrat;

        /* correct for pressure difference from standard pressure */
        tsol[iw] = tsol[iw]*exp(-0.5*taur[iw]/mu0*(pr/p0-1));
        tsen[iw] = tsen[iw]*exp(-0.5*taur[iw]/mu *(pr/p0-1));

        if ((evalmask & TRANSSPHER) != 0) {
            /* correct for airmass difference, plane-parallel to spherical atmosphere */
            tsol[iw] = pow(tsol[iw],airmass_sph/airmass_plp);
            tsen[iw] = pow(tsen[iw],airmass_sph/airmass_plp);
        }
    }

    return;
}


/* ---------------------------------------------------------------------------------------- */
/* wangaer() - compute aerosol reflectance using Gordon & Wang 1994 algorithm               */
/*                                                                                          */
/* B. Franz, 1 June 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
int wangaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
        int32_t nmodels, int32_t mindx[],
        float solz, float senz, float phi, float wv, float rhoa[],
        int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir)
{
    int    modflg;
    float *epsmin;
    float *epsmax;

    float epsmin1;
    float epsmax1;

    float rhoasmin[NBANDS];
    float rhoasmax[NBANDS];
    float rhoamin [NBANDS];
    float rhoamax [NBANDS];

    float cc = 0.0;
    int   iw;


    /* find upper and lower-bounding models */
    modflg = model_select_wang2(sensorID,wave,nwave,nmodels,mindx,solz,senz,phi,wv,rhoa,iwnir_s,iwnir_l,
            modmin,modmax,modrat,epsnir);

    //    modflg = model_select_wang(sensorID,wave,nwave,solz,senz,phi,wv,rhoa,iwnir_s,iwnir_l,
    //                               modmin,modmax,modrat,epsnir);

    /* if no lower-bounding aerosol model, set-up for extrapolation */
    if (modflg < 0)
        cc = log(*epsnir)/(wave[iwnir_l]-wave[iwnir_s]);

    /* get model epsilon for each bounding model, all wavelengths */
    epsmin = model_epsilon(*modmin,iwnir_l,wave,nwave,solz,senz,phi);
    epsmax = model_epsilon(*modmax,iwnir_l,wave,nwave,solz,senz,phi);

    /* get SS aerosol reflectance at longest wavelength for the two models */
    if (rhoa_to_rhoas(sensorID,*modmin,solz,senz,phi,wv,rhoa,wave,nwave,iwnir_l,iwnir_l,rhoasmin) != 0)
        return(1);
    if (rhoa_to_rhoas(sensorID,*modmax,solz,senz,phi,wv,rhoa,wave,nwave,iwnir_l,iwnir_l,rhoasmax) != 0)
        return(1);

    /* compute SS aerosol reflectance in all bands */
    for (iw=0; iw<nwave; iw++) {

        epsmin1 = epsmin[iw];
        epsmax1 = epsmax[iw];

        if (modflg < 0) {
            epsmin1 = exp(cc*(wave[iwnir_l]-wave[iw]));
            epsmax1 = epsmin1;
        }

        rhoasmin[iw] = rhoasmin[iwnir_l]*epsmin1;
        rhoasmax[iw] = rhoasmax[iwnir_l]*epsmax1;
    }

    /* compute MS aerosol reflectance in visible bands */
    rhoas_to_rhoa(sensorID,*modmin,solz,senz,phi,wv,rhoasmin,wave,nwave,0,nwave-1,rhoamin);
    rhoas_to_rhoa(sensorID,*modmax,solz,senz,phi,wv,rhoasmax,wave,nwave,0,nwave-1,rhoamax);

    /* interpolate between upper and lower-bounding models */
    for (iw=0; iw<nwave; iw++) {
        rhoa[iw] = rhoamin[iw]*(1.0-(*modrat)) + rhoamax[iw]*(*modrat);
    }

    return(0);
}

/* ---------------------------------------------------------------------------------------- */
/* model_select_franz() - Franz aerosol model selection process.                            */
/*                                                                                          */
/* B. Franz, 1 February 2009.                                                               */
/* ---------------------------------------------------------------------------------------- */

typedef struct rhoaT_struct {
    int32_t   modnum;
    float  rhoa;
    float  eps;
} rhoaTstr;

int comp_rhoaT(rhoaTstr *x, rhoaTstr *y) {return(x->rhoa < y->rhoa ? -1 : 1);}

int model_select_franz(int32_t sensorID, float wave[], int32_t nwave, int32_t nmodel, int32_t mindx[],
        float solz, float senz, float phi,
        float wv, float rhoa[], int32_t iwnir_s, int32_t iwnir_l,
        int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir)
{

    float rhoas[NBANDS];
    float rhoa_tmp[NBANDS];
    rhoaTstr rhoa_tab[MAXMODEL];

    float *eps;
    int   jm, im;
    int   jm1, jm2;
    float wt;

    *modmin =  -1;
    *modmax =  -1;
    *modrat = 0.0;
    *epsnir = 0.0;

    // predict MS aerosol reflectance assuming each model

    for (jm=0; jm<nmodel; jm++) {

        im = mindx[jm];

        // get model epsilon at this geometry
        eps = model_epsilon(im,iwnir_l,wave,nwave,solz,senz,phi);

        // get SS aerosol reflectance at iwnir_l
        if (rhoa_to_rhoas(sensorID,im,solz,senz,phi,wv,rhoa,wave,nwave,iwnir_l,iwnir_l,rhoas) != 0)
            return(1);

        // get SS aerosol reflectance at iwnir_s
        rhoas[iwnir_s] = rhoas[iwnir_l] * eps[iwnir_s];

        // get MS aerosol reflectance at iwnir_l and save
        rhoas_to_rhoa(sensorID,im,solz,senz,phi,wv,rhoas,wave,nwave,iwnir_s,iwnir_s,rhoa_tmp);

        rhoa_tab[jm].modnum = im;
        rhoa_tab[jm].rhoa = rhoa_tmp[iwnir_s];
        rhoa_tab[jm].eps  = eps[iwnir_s];
    }

    // put results in ascending order of predicted rhoa[iwnir_s]
    qsort(rhoa_tab,nmodel,sizeof(rhoaTstr),(int (*)(const void *,const void *)) comp_rhoaT);

    // compare observed rhoa with model predictions at iwnir_s to select models
    for (jm=0; jm<nmodel; jm++) {
        if (rhoa_tab[jm].rhoa > rhoa[iwnir_s])
            break;
    }
    if (jm == 0) {
        jm1 = 0;
        jm2 = 1;
    } else if (jm == nmodel) {
        jm1 = nmodel-2;
        jm2 = nmodel-1;
    } else {
        jm1 = jm  - 1;
        jm2 = jm1 + 1;
    }
    wt  = (rhoa[iwnir_s] - rhoa_tab[jm1].rhoa)/(rhoa_tab[jm2].rhoa - rhoa_tab[jm1].rhoa);

    *modmin = rhoa_tab[jm1].modnum;
    *modmax = rhoa_tab[jm2].modnum;
    *modrat = wt;
    *epsnir = rhoa_tab[jm1].eps*(1.0-wt) + rhoa_tab[jm2].eps*wt;

    return(0);
}

/* ---------------------------------------------------------------------------------------- */
/* franzaer() - compute aerosol reflectance using modified G&W algorithm                    */
/*                                                                                          */
/* B. Franz, February 2009.                                                                 */
/* ---------------------------------------------------------------------------------------- */
int franzaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
        int32_t nmodels, int32_t mindx[],
        float solz, float senz, float phi, float wv, float rhoa[],
        int32_t *modmin, int32_t *modmax, float *modrat, float *epsnir)
{
    float *epsmin;
    float *epsmax;

    float epsmin1;
    float epsmax1;

    float rhoasmin[NBANDS];
    float rhoasmax[NBANDS];
    float rhoamin [NBANDS];
    float rhoamax [NBANDS];

    int   iw;


    /* find upper and lower-bounding models */
    if (model_select_franz(sensorID,wave,nwave,nmodels,mindx,solz,senz,phi,wv,rhoa,iwnir_s,iwnir_l,
                modmin,modmax,modrat,epsnir) != 0)
        return(1);

    /* get model epsilon for each bounding model, all wavelengths */
    epsmin = model_epsilon(*modmin,iwnir_l,wave,nwave,solz,senz,phi);
    epsmax = model_epsilon(*modmax,iwnir_l,wave,nwave,solz,senz,phi);

    /* get SS aerosol reflectance at longest wavelength for the two models */
    if (rhoa_to_rhoas(sensorID,*modmin,solz,senz,phi,wv,rhoa,wave,nwave,iwnir_l,iwnir_l,rhoasmin) != 0)
        return(1);
    if (rhoa_to_rhoas(sensorID,*modmax,solz,senz,phi,wv,rhoa,wave,nwave,iwnir_l,iwnir_l,rhoasmax) != 0)
        return(1);

    /* compute SS aerosol reflectance in all bands */
    for (iw=0; iw<nwave; iw++) {

        epsmin1 = epsmin[iw];
        epsmax1 = epsmax[iw];

        rhoasmin[iw] = rhoasmin[iwnir_l]*epsmin1;
        rhoasmax[iw] = rhoasmax[iwnir_l]*epsmax1;
    }

    /* compute MS aerosol reflectance in visible bands */
    rhoas_to_rhoa(sensorID,*modmin,solz,senz,phi,wv,rhoasmin,wave,nwave,0,nwave-1,rhoamin);
    rhoas_to_rhoa(sensorID,*modmax,solz,senz,phi,wv,rhoasmax,wave,nwave,0,nwave-1,rhoamax);

    /* interpolate between upper and lower-bounding models */
    for (iw=0; iw<nwave; iw++) {
        rhoa[iw] = rhoamin[iw]*(1.0-(*modrat)) + rhoamax[iw]*(*modrat);
    }

    return(0);
}

/* ---------------------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------------------------- */
static int order_models(const void *p1, const void *p2)
{
    aermodstr *x = *(aermodstr **)p1;
    aermodstr *y = *(aermodstr **)p2;

    if (x->rh == y->rh) {
        if (x->sd > y->sd)
            return( 1);
        else
            return(-1);
    } else {
        if (x->rh > y->rh)
            return( 1);
        else
            return(-1);
    }
}


/* ---------------------------------------------------------------------------------------- */
/* rhaer() - compute aerosol reflectance using RH selection + G & W 1994                    */
/*                                                                                          */
/* B. Franz, 1 June 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
int rhaer(int32_t sensorID, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
        float solz, float senz, float phi, float wv, float rh, float pr, float taur[], float rhoa[],
        int32_t *modmin, int32_t *modmax, float *modrat, float *eps, float taua[], float tsol[], float tsen[])
{
    static int   firstCall = 1;
    static int   nrh;
    static float rhtab[MAXAERMOD];
    static int   nsd;
    static int   sdtab[MAXAERMOD];

    float rhoa1[NBANDS];
    float rhoa2[NBANDS];
    float taua1[NBANDS];
    float taua2[NBANDS];
    float tsol1[NBANDS];
    float tsol2[NBANDS];
    float tsen1[NBANDS];
    float tsen2[NBANDS];
    float eps1;
    float eps2;
    float modrat1;
    float modrat2;
    int32_t  modmin1;
    int32_t  modmin2;
    int32_t  modmax1;
    int32_t  modmax2;
    int32_t  mindx1[MAXAERMOD];
    int32_t  mindx2[MAXAERMOD];
    int   irh1, irh2, irh;
    int   isd;
    float wt;

    int   iw, im;

    if (firstCall) {
        firstCall = 0;
        float lastrh = -1.0;
        int   lastsd = -1;
        if (!have_rh || !have_sd) {
            HY1C_out("-E- %s line %d: This aerosol selection method requires models with a Relative Humidity attribute.\n",
                    __FILE__,__LINE__);
            exit(1);
        }
        // need in order of rh and sd within rh
        qsort(aertab->model,aertab->nmodel,sizeof(aermodstr*),(int (*)(const void *,const void *)) order_models);
        nsd = 0;
        nrh = 0;
        for (im=0; im<aertab->nmodel; im++) {
            if (aertab->model[im]->rh != lastrh) {
                rhtab[nrh] = aertab->model[im]->rh;
                lastrh     = rhtab[nrh];
                nrh++;
            }
            if (nrh == 1 && aertab->model[im]->sd != lastsd) {
                sdtab[nsd] = aertab->model[im]->sd;
                lastsd     = sdtab[nsd];
                nsd++;
            }
        }
        if (nrh*nsd != aertab->nmodel) {
            HY1C_out("-E- %s line %d: number of humidities (%d) x number of size distributions (%d) must equal number of models (%d).\n",
                    __FILE__,__LINE__,nrh,nsd,aertab->nmodel);
            exit(1);
        } else {
            HY1C_out("%d aerosol models: %d humidities x %d size fractions\n",aertab->nmodel,nrh,nsd);
            for (irh=0; irh<nrh; irh++) {
                for (isd=0; isd<nsd; isd++) {
                    im = irh*nsd+isd;
                    HY1C_out("model %d, rh=%f, sd=%d, alpha=%f, name=%s\n",
                            im,aertab->model[im]->rh, aertab->model[im]->sd, aertab->model[im]->angstrom[1],aertab->model[im]->name);
                }
            }
        }
    }

    // initialize

    *modmin = -1;
    *modmax = -1;
    *modrat =  0;

    for (iw=0; iw<nwave; iw++) {
        taua[iw] = -1.0;
        tsol[iw] = -1.0;
        tsen[iw] = -1.0;
        rhoa1[iw] = rhoa[iw];
        rhoa2[iw] = rhoa[iw];
        rhoa [iw] = BAD_FLT;
    }

    // find RH index and wts

    if (nrh == 1 || rhtab[0] > rh) {
        irh1 = 0;
        irh2 = 0;
        wt = 0.0;
    } else if (rhtab[nrh-1] < rh) {
        irh1 = nrh-1;
        irh2 = nrh-1;
        wt = 0.0;
    } else {
        for (irh=0; irh<nrh; irh++) {
            if (rhtab[irh] > rh)
                break;
        }
        irh1 = MIN(MAX(0,irh-1),nrh-2);
        irh2 = irh1 + 1;
        wt = (rh - rhtab[irh1])/(rhtab[irh2]-rhtab[irh1]);
    }

    // set indices of active model sets

    for (im=0; im<nsd; im++) {
        mindx1[im] = irh1*nsd + im;
        mindx2[im] = irh2*nsd + im;
    }

    // find compute aerosol reflectances, aot, diffuse trans, eps from first model set

    if (aer_opt == AERRHFRNIR) {
        if ( franzaer(sensorID,wave,nwave,iwnir_s,iwnir_l,nsd,mindx1,
                    solz,senz,phi,wv,rhoa1,&modmin1,&modmax1,&modrat1,&eps1) != 0) {
            return(1);
        }
    } else {
        if ( wangaer(sensorID,wave,nwave,iwnir_s,iwnir_l,nsd,mindx1,
                    solz,senz,phi,wv,rhoa1,&modmin1,&modmax1,&modrat1,&eps1) != 0) {
            return(1);
        }
    }

    diff_tran(sensorID,wave,nwave,iwnir_l,solz,senz,phi,wv,pr,taur,
            modmin1,modmax1,modrat1,rhoa1,taua1,tsol1,tsen1);

    // compute aerosol reflectances, aot, diffuse trans, eps from second model set (if needed)

    if (irh2 != irh1) {

        if (aer_opt == AERRHFRNIR) {
            if ( franzaer(sensorID,wave,nwave,iwnir_s,iwnir_l,nsd,mindx2,
                        solz,senz,phi,wv,rhoa2,&modmin2,&modmax2,&modrat2,&eps2) != 0) {
                return(1);
            }
        } else {
            if ( wangaer(sensorID,wave,nwave,iwnir_s,iwnir_l,nsd,mindx2,
                        solz,senz,phi,wv,rhoa2,&modmin2,&modmax2,&modrat2,&eps2) != 0) {
                return(1);
            }
        }

        diff_tran(sensorID,wave,nwave,iwnir_l,solz,senz,phi,wv,pr,taur,
                modmin2,modmax2,modrat2,rhoa2,taua2,tsol2,tsen2);

        for (iw=0; iw<nwave; iw++) {
            rhoa[iw] = rhoa1[iw]*(1-wt) + rhoa2[iw]*wt;
            taua[iw] = taua1[iw]*(1-wt) + taua2[iw]*wt;
            tsol[iw] = tsol1[iw]*(1-wt) + tsol2[iw]*wt;
            tsen[iw] = tsen1[iw]*(1-wt) + tsen2[iw]*wt;
        }
        *eps = eps1*(1-wt) + eps2*wt;

    } else {

        for (iw=0; iw<nwave; iw++) {
            rhoa[iw] = rhoa1[iw];
            taua[iw] = taua1[iw];
            tsol[iw] = tsol1[iw];
            tsen[iw] = tsen1[iw];
        }
        *eps = eps1;
    }

    // store model info for the dominat RH only

    if (wt < 0.5) {
        *modmin = modmin1;
        *modmax = modmax1;
        *modrat = modrat1;
    } else {
        *modmin = modmin2;
        *modmax = modmax2;
        *modrat = modrat2;
    }

    return(0);
}


/* ---------------------------------------------------------------------------------------- */
/* fixedaer() - compute aerosol reflectance for fixed aerosol model                         */
/*                                                                                          */
/* B. Franz, August 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
int fixedaer(int32_t sensorID, int32_t modnum, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
        char models[MAXAERMOD][32], int32_t nmodels,
        float solz, float senz, float phi, float wv, float rhoa[],float *epsnir)
{
    float *eps;
    float  rhoas[NBANDS];
    int    iw;

    if (rhoa[iwnir_l] < 0.0) {
        *epsnir = BAD_FLT;
        for (iw=0; iw<nwave; iw++)
            rhoas[iw] = BAD_FLT;
        return(1);
    }

    /* get model epsilon for all wavelengths at this geometry */
    eps = model_epsilon(modnum,iwnir_l,wave,nwave,solz,senz,phi);

    /* get SS aerosol reflectance at longest wavelength */
    if (rhoa_to_rhoas(sensorID,modnum,solz,senz,phi,wv,rhoa,wave,nwave,iwnir_l,iwnir_l,rhoas) != 0) {
        HY1C_out("Error getting rhoas\n");
        return(1);
    }

    /* compute SS aerosol reflectance in visible bands */
    for (iw=0; iw<nwave; iw++) {
        rhoas[iw] = rhoas[iwnir_l]*eps[iw];
    }

    /* compute MS aerosol reflectance in visible bands */
    rhoas_to_rhoa(sensorID,modnum,solz,senz,phi,wv,rhoas,wave,nwave,0,nwave-1,rhoa);

    if (iwnir_s == iwnir_l)
        *epsnir = eps[iwnir_l-1];
    else
        *epsnir = eps[iwnir_s];

    return(0);
}


/* ---------------------------------------------------------------------------------------- */
/* fixedmodpair() - compute aerosol reflectance for fixed model pair                        */
/*                                                                                          */
/* B. Franz, August 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
int fixedmodpair(
        int32_t sensorID, instr *input, float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
        float solz, float senz, float phi, float wv, int32_t modmin, int32_t modmax, float modrat,
        float rhoa[],float *eps)
{
    float rhoa1 [NBANDS];
    float rhoa2 [NBANDS];
    float eps1;
    float eps2;
    int   iw;
    int   status;

    if (modmin < 0 || modmin >= input->naermodels ||
            modmax < 0 || modmax >= input->naermodels ||
            modrat < 0.0 || modrat > 1.0) {
        HY1C_out("Bogus input for fixed model pair: %d %d %f\n",modmin+1,modmax+1,modrat);
        exit(1);
    }


    if  (rhoa[iwnir_l] > input->rhoamin) {

        rhoa2[iwnir_l] = rhoa1[iwnir_l] = rhoa[iwnir_l];

        /* get aerosol reflectance in visible for fixed model, and epsilon */
        status = fixedaer(sensorID,modmin,wave,nwave,iwnir_s,iwnir_l,
                input->aermodels,input->naermodels,
                solz,senz,phi,wv,rhoa1,&eps1);
        status = fixedaer(sensorID,modmax,wave,nwave,iwnir_s,iwnir_l,
                input->aermodels,input->naermodels,
                solz,senz,phi,wv,rhoa2,&eps2);

        /* convert aerosol relectance to radiance */
        if (status == 0) {
            for (iw=0; iw<nwave; iw++) {
                if (iw != iwnir_l) // without this check tLw-La may go slight negative
                    rhoa[iw] = MAX((1.0-modrat)*rhoa1[iw] + modrat*rhoa2[iw],0.0);
            }
            *eps  = (1.0-modrat)*eps1 + modrat*eps2;
        }

    } else if (rhoa[iwnir_l] > -(input->rhoamin)) {

        /* if input NIR is near zero, assume a small white aerosol signal */
        *eps    = 1.0;
        for (iw=0; iw<nwave; iw++) {
            rhoa[iw] = MAX(rhoa[iwnir_l],1e-6);
        }

        status = 0;

    } else {

        /* if input NIR is very negative, fail the case */
        status = 1;
    }


    return(status);
}


/* ---------------------------------------------------------------------------------------- */
/* fixedaot() - compute aerosol reflectance for fixed aot(lambda)                           */
/*                                                                                          */
/* B. Franz, August 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
int fixedaot(int32_t sensorID, float aot[], float wave[], int32_t nwave, int32_t iwnir_s, int32_t iwnir_l,
        float solz, float senz, float phi, float wv, int32_t *modmin, int32_t *modmax, float *modrat,
        float rhoa[],float *epsnir)
{
    static int firstCall = 1;
    static int angst_band1 = -1;
    static int angst_band2 = -1;

    float *phase1;
    float *phase2;
    float f1    [MAXBAND];
    float f2    [MAXBAND];
    float lnf1  [MAXBAND];
    float lnf2  [MAXBAND];
    float rhoas1[NBANDS];
    float rhoas2[NBANDS];
    float rhoa1 [NBANDS];
    float rhoa2 [NBANDS];
    float eps1;
    float eps2;
    float angstrom;
    int   iw, iwtab;

    if (firstCall) {
        angst_band1 = windex(520,wave,nwave);
        angst_band2 = windex(865,wave,nwave);
        firstCall = 0;
    }

    /* bail on negative input aot */
    for (iw=0; iw<nwave; iw++) {
        if (aot[iw] < 0.0) return(1);
    }

    /* compute angstrom and use to select bounding models */
    if (aot[iwnir_l] > 0.0)
        angstrom = -log(aot[angst_band1]/aot[angst_band2])/
            log(wave[angst_band1]/wave[angst_band2]);
    else
        angstrom = 0.0;

    model_select_angstrom(angstrom,modmin,modmax,modrat);


    /* get model phase function for all wavelengths at this geometry for the two models */
    phase1 = model_phase(*modmin,solz,senz,phi);
    phase2 = model_phase(*modmax,solz,senz,phi);

    /* compute factor for SS approximation, set-up for interpolation */
    for (iw=0; iw<aertab->nwave; iw++) {
        f1[iw] = aertab->model[*modmin]->albedo[iw]*phase1[iw]/4.0/mu0/mu;
        f2[iw] = aertab->model[*modmax]->albedo[iw]*phase2[iw]/4.0/mu0/mu;
        if (interp) {
            lnf1[iw] = log(f1[iw]);
            lnf2[iw] = log(f2[iw]);
        }
    }

    /* compute SS aerosol reflectance */
    if (interp) {
        for (iw=0; iw<nwave; iw++) {
            iwtab = iwatab[iw];
            if (aertab->wave[iwtab] != wave[iw] && wave[iw] > 0) {
                rhoas1[iw] = aot[iw]*exp(linterp(aertab->wave,lnf1,aertab->nwave,wave[iw]));
                rhoas2[iw] = aot[iw]*exp(linterp(aertab->wave,lnf2,aertab->nwave,wave[iw]));
            } else {
                rhoas1[iw] = aot[iw]*f1[iwtab];
                rhoas2[iw] = aot[iw]*f2[iwtab];
            }
        }
    } else {
        for (iw=0; iw<nwave; iw++) {
            iwtab = iwatab[iw];
            rhoas1[iw] = aot[iw]*f1[iwtab];
            rhoas2[iw] = aot[iw]*f2[iwtab];
        }
    }
    eps1 = rhoas1[iwnir_s]/rhoas1[iwnir_l];
    eps2 = rhoas2[iwnir_s]/rhoas2[iwnir_l];

    /* compute MS aerosol reflectance */
    rhoas_to_rhoa(sensorID,*modmin,solz,senz,phi,wv,rhoas1,wave,nwave,0,nwave-1,rhoa1);
    rhoas_to_rhoa(sensorID,*modmax,solz,senz,phi,wv,rhoas2,wave,nwave,0,nwave-1,rhoa2);

    for (iw=0; iw<nwave; iw++) {
        rhoa[iw] = (1.0-*modrat) * rhoa1[iw] + *modrat * rhoa2[iw];
    }
    *epsnir = (1.0-*modrat) * eps1 + *modrat * eps2;

    return(0);
}


/* ---------------------------------------------------------------------------------------- */
/* aerosol() - compute aerosol reflectance using specified algorithm                        */
/*                                                                                          */
/* B. Franz, 1 June 2004.                                                                   */
/* ---------------------------------------------------------------------------------------- */
int aerosol(l2str *l2rec, int32_t aer_opt_in, aestr *aerec, int32_t ip,
        float wave[], int32_t nwave, int32_t iwnir_s_in, int32_t iwnir_l_in,
        float F0_in[], float La1_in[], float La2_out[],
        float t_sol_out[], float t_sen_out[], float *eps, float taua_out[],
        int32_t *modmin, int32_t *modmax, float *modrat)
{
    static int firstCall = 1;
    static int32_t mindx[MAXAERMOD];

    int   status = 1;
    float rhoa   [NBANDS];
    float radref [NBANDS];
    float temp;
    float *aot;
    float angstrom;

    int iw;

    float F0   [NBANDS];
    float taur [NBANDS];
    float La1  [NBANDS];
    float La2  [NBANDS];
    float t_sol[NBANDS];
    float t_sen[NBANDS];
    float taua [NBANDS];

    int32_t  sensorID = l2rec->sensorID;
    float solz     = l2rec->solz  [ip];
    float senz     = l2rec->senz  [ip];
    float phi      = l2rec->delphi[ip];
    float wv       = l2rec->wv    [ip];
    float rh       = l2rec->rh    [ip];
    float pr       = l2rec->pr    [ip];

    instr *input = l2rec->input;

    /* set static global evaluation level */
    evalmask = input->evalmask;
    aer_opt  = aer_opt_in;

    /* transfer inputs per band to inputs per wavelength */
    for (iw=0; iw<nwave; iw++) {
        F0  [iw] = F0_in [iw];
        taur[iw] = l2rec->Tau_r[iw];
        La1 [iw] = La1_in[iw];
        if (iw == iwnir_s_in) iwnir_s = iw;
        if (iw == iwnir_l_in) iwnir_l = iw;
    }

    /* compute total airmass (static global) */
    mu0     = cos(solz/radeg);
    mu      = cos(senz/radeg);
    airmass = 1.0/mu0 + 1.0/mu;
    if ((evalmask & TRANSSPHER) != 0) {
        airmass_plp = pp_airmass(solz) + pp_airmass(senz);
        airmass_sph = ky_airmass(solz) + ky_airmass(senz);
    }

    /* initialize epsilon and diffuse transmittance to defaults */
    *eps = 1.0;
    *modmin = -1;
    *modmax = -1;
    *modrat = -1;
    for (iw=0; iw<nwave; iw++) {
        taua  [iw] = 0.0;
        t_sol [iw] = 1.0;
        t_sen [iw] = 1.0;
        La2   [iw] = 0.0;
        radref[iw] = pi / F0[iw] / mu0;
    }

    /* load aerosol model tables */
    if (!loaded) {
        int32_t im;
        load_aermod(sensorID,wave,nwave,input->aermodels,input->naermodels, NULL);
        for (im=0; im<aertab->nmodel; im++) mindx[im] = im;
        if (have_rh && aertab->nmodel >= 30) {
            HY1C_out("\nLimiting aerosol models based on RH.\n");
            use_rh = 1;
        }
    }

    /* Change the aerosol option if need be */
    if (use_rh)
        switch (aer_opt) {
            case AERWANG:     aer_opt = AERRH;  break;
            case AERWANGNIR:  aer_opt = AERRHNIR;  break;
            case AERWANGSWIR: aer_opt = AERRHSWIR; break;
            case AERMUMM:     aer_opt = AERRHMUMM; break;
        }


    if (firstCall) {
        if (aer_opt > 0 && aer_opt <= MAXAERMOD) {
            HY1C_out("\nUsing fixed aerosol model #%d (%s)\n",aer_opt,input->aermodels[aer_opt-1]);
            HY1C_out("Extrapolating from %4.1f nm\n",wave[iwnir_l]);
        } else {
            switch (aer_opt) {
                case AERWHITE:
                    HY1C_out("\nUsing white-aerosol approximation\n");
                    HY1C_out("Extrapolating from %4.1f nm\n",wave[iwnir_l]);
                    break;
                case AERWANG:
                case AERRH:
                    HY1C_out("\nUsing Gordon & Wang aerosol model selection\n");
                    HY1C_out("Using bands at %4.1f and %4.1f nm for model selection\n",wave[iwnir_s],wave[iwnir_l]);
                    HY1C_out("Extrapolating from %4.1f nm\n",wave[iwnir_l]);
                    break;
                case AERWANGNIR:
                case AERRHNIR:
                    HY1C_out("\nUsing Gordon & Wang aerosol model selection\n");
                    HY1C_out("  and NIR correction with up to %d iterations\n",input->aer_iter_max);
                    HY1C_out("Using bands at %4.1f and %4.1f nm for model selection\n",wave[iwnir_s],wave[iwnir_l]);
                    HY1C_out("Extrapolating from %4.1f nm\n",wave[iwnir_l]);
                    break;
                case AERWANGSWIR:
                case AERRHSWIR:
                    HY1C_out("\nUsing Gordon & Wang aerosol model selection with NIR/SWIR switching (only SWIR).\n");
                    /*HY1C_out("NIR correction with up to %d iterations\n",input->aer_iter_max);
                      HY1C_out("NIR bands at %d and %d nm\n",input->aer_wave_short,input->aer_wave_long);*/
                    HY1C_out("SWIR bands at %d and %d nm\n\n",input->aer_swir_short, input->aer_swir_long);
                    break;
                case AERMUMM:
                case AERRHMUMM:
                    HY1C_out("\nUsing Gordon & Wang aerosol model selection\n");
                    HY1C_out("  and MUMM correction\n");
                    HY1C_out("Using bands at %4.1f and %4.1f nm for model selection\n",wave[iwnir_s],wave[iwnir_l]);
                    HY1C_out("Extrapolating from %4.1f nm\n",wave[iwnir_l]);
                    break;
                case AERRHFRNIR:
                    HY1C_out("\nUsing Ahmad & Franz aerosol model selection.\n");
                    HY1C_out("  and NIR correction with up to %d iterations\n",input->aer_iter_max);
                    HY1C_out("Using bands at %4.1f and %4.1f nm for model selection\n",wave[iwnir_s],wave[iwnir_l]);
                    HY1C_out("Extrapolating from %4.1f nm\n",wave[iwnir_l]);
                    break;
                case FIXAOT:
                    HY1C_out("\nUsing fixed, input aerosol optical thickenesses for aerosol selection.\n");
                    break;
                case FIXMODPAIR:
                    HY1C_out("\nUsing multiple scattering aerosols with fixed model pair\n");
                    break;
                case FIXMODPAIRNIR:
                    HY1C_out("\nUsing multiple scattering aerosols with fixed model pair\n");
                    HY1C_out("  and NIR iteration with up to %d iterations\n",input->aer_iter_max);
                    HY1C_out("Extrapolating from %4.1f nm\n",wave[iwnir_l]);
                    break;
                case FIXANGSTROM:
                    HY1C_out("\nUsing fixed aerosol model based on predefined Angstrom exponent\n");
                    HY1C_out("Extrapolating from %4.1f nm\n",wave[iwnir_l]);
                    break;
                case FIXANGSTROMNIR:
                    HY1C_out("\nUsing fixed aerosol model based on predefined Angstrom exponent\n");
                    HY1C_out("  and NIR iteration with up to %d iterations\n",input->aer_iter_max);
                    HY1C_out("Extrapolating from %4.1f nm\n",wave[iwnir_l]);
                    break;
                default:
                    HY1C_out("\nErroneous aerosol option, %d\n",aer_opt);
                    exit(FATAL_ERROR);
                    break;
            }
        }
        firstCall = 0;
    }


    /* convert input aerosol radiances to relectance */
    for (iw=iwnir_s; iw<=iwnir_l; iw+=MAX(iwnir_l-iwnir_s,1))
        rhoa[iw] = La1[iw] * radref[iw];


    /* compute aerosol using requested method */
    /* -------------------------------------- */

    switch (aer_opt) {

        case AERWANG: case AERWANGNIR: case AERWANGSWIR: case AERMUMM:

            if (iwnir_l <= iwnir_s || wave[iwnir_s] < 600) {
                HY1C_out("Aerosol selection bands must be greater than 600nm with short wave less than long wave (%d,%d)\n",iwnir_l,iwnir_s);
                exit(1);
            }

            /* Multi-Scattering with Gordon & Wang Aerosol Selection */
            /* ----------------------------------------------------- */

            /* convert input NIR aerosol radiances to relectance */
            for (iw=iwnir_s; iw<=iwnir_l; iw++) {
                rhoa[iw] = La1[iw] * radref[iw];
            }

            /* require sufficient signal in two NIR bands */
            if  (rhoa[iwnir_s] > input->rhoamin && rhoa[iwnir_l] > input->rhoamin) {

                /* require MS epsilon to be reasonable */
                if (La1[iwnir_s]/La1[iwnir_l] > 0.1) {

                    status = wangaer(sensorID,wave,nwave,iwnir_s,iwnir_l,
                            aertab->nmodel,mindx,
                            solz,senz,phi,wv,rhoa,modmin,modmax,modrat,eps);

                    if (status == 0)
                        for (iw=0; iw<nwave; iw++)
                            La2[iw] = rhoa[iw] / radref[iw];
                }

            } else if (rhoa[iwnir_s] > -(input->rhoamin) && rhoa[iwnir_l] > -(input->rhoamin)) {

                /* if input NIR is near zero, assume a small white aerosol signal */
                *eps    = 1.0;
                *modmin = 0;
                *modmax = 0;
                *modrat = 0.0;
                temp    = MAX(rhoa[iwnir_l],1e-6);
                for (iw=0; iw<nwave; iw++) {
                    rhoa[iw] = temp;
                    La2 [iw] = rhoa[iw] / radref[iw];
                }

                status = 0;

            } else {

                /* if input NIR is very negative, fail the case */
                status = 1;
            }

            break;

        case AERRH:
        case AERRHNIR:
        case AERRHSWIR:
        case AERRHMUMM:
        case AERRHFRNIR:

            if (iwnir_l <= iwnir_s || wave[iwnir_s] < 600) {
                HY1C_out("Aerosol selection bands must be greater than 600nm with short wave less than long wave");
                exit(1);
            }

            /* Multi-Scattering with RH-based selection              */
            /* ----------------------------------------------------- */

            /* convert input NIR aerosol radiances to relectance */
            for (iw=iwnir_s; iw<=iwnir_l; iw++) {
                rhoa[iw] = La1[iw] * radref[iw];
            }

            /* require sufficient signal in two NIR bands */
            if  (rhoa[iwnir_s] > input->rhoamin && rhoa[iwnir_l] > input->rhoamin) {

                /* require MS epsilon to be reasonable */
                if (rhoa[iwnir_s]/rhoa[iwnir_l] > 0.1 && rhoa[iwnir_s]/rhoa[iwnir_l] < 10.0) {

                    status = rhaer(sensorID,wave,nwave,iwnir_s,iwnir_l,
                            solz,senz,phi,wv,rh,pr,taur,rhoa,
                            modmin,modmax,modrat,eps,taua,t_sol,t_sen);

                    if (status == 0)
                        for (iw=0; iw<nwave; iw++)
                            La2[iw] = rhoa[iw] / radref[iw];
                }

            } else if (rhoa[iwnir_s] > -(input->rhoamin) && rhoa[iwnir_l] > -(input->rhoamin)) {

                /* if input NIR is near zero, assume a small white aerosol signal */
                *eps    = 1.0;
                *modmin = aertab->nmodel;
                *modmax = aertab->nmodel;
                *modrat = 0.0;
                temp    = MAX(rhoa[iwnir_l],1e-6);
                for (iw=0; iw<nwave; iw++) {
                    rhoa[iw] = temp;
                    La2 [iw] = rhoa[iw] / radref[iw];
                }

                diff_tran(sensorID,wave,nwave,iwnir_l,solz,senz,phi,wv,pr,taur,
                        *modmin,*modmax,*modrat,rhoa,taua,t_sol,t_sen);

                status = 0;

            } else {

                /* if input NIR is very negative, fail the case */
                status = 1;
            }

            break;

        case AERWHITE:

            /* White Aerosol */
            /* ------------- */

            if ( La1[iwnir_l] > 0.0) {

                *eps    = 1.0;
                *modmin = 0;
                *modmax = 0;
                *modrat = 0.0;

                for (iw=0; iw<nwave; iw++) {
                    La2 [iw] = *eps * F0[iw]/F0[iwnir_l]*La1[iwnir_l];
                    rhoa[iw] = La2[iw] * radref[iw];
                }


                status = 0;
            }
            break;

        case FIXMODPAIR: case FIXMODPAIRNIR:

            /* Multi-Scattering with Fixed Model Pair */
            /* -------------------------------------- */

            if (aerec != NULL && aerec->mode == ON) {
                *modmin = aerec->mod_min[ip] - 1;
                *modmax = aerec->mod_max[ip] - 1;
                *modrat = aerec->mod_rat[ip];
            } else {
                *modmin = input->aermodmin - 1;
                *modmax = input->aermodmax - 1;
                *modrat = input->aermodrat;
            }

            status = fixedmodpair(sensorID,input,wave,nwave,iwnir_s,iwnir_l,solz,senz,phi,wv,
                    *modmin,*modmax,*modrat,rhoa,eps);
            if (status == 0) {
                for (iw=0; iw<nwave; iw++)
                    La2[iw] = rhoa[iw] / radref[iw];
            }

            break;

        case FIXANGSTROM: case FIXANGSTROMNIR:

            if (input->aer_angstrom > -2) {
                angstrom = input->aer_angstrom;
            } else {
                angstrom = bin_climatology(input->aerbinfile,*(l2rec->day),l2rec->lon[ip],l2rec->lat[ip],"angstrom");
            }

            if (angstrom > -2) {

                model_select_angstrom(angstrom,modmin,modmax,modrat);

                status = fixedmodpair(sensorID,input,wave,nwave,iwnir_s,iwnir_l,solz,senz,phi,wv,
                        *modmin,*modmax,*modrat,rhoa,eps);
            } else
                status = 1;

            if (status == 0) {
                for (iw=0; iw<nwave; iw++)
                    La2[iw] = rhoa[iw] / radref[iw];
            }

            break;

        case FIXAOT:

            /* Multi-Scattering with fixed AOT */
            /* ------------------------------- */
            if (aerec != NULL && aerec->mode == ON)
                aot = &aerec->taua[ip*NBANDS];
            else
                aot = input->taua;

            status = fixedaot(sensorID,aot,wave,nwave,iwnir_s,iwnir_l,solz,senz,phi,wv,
                    modmin,modmax,modrat,rhoa,eps);
            if (status == 0) {
                for (iw=0; iw<nwave; iw++)
                    La2[iw] = rhoa[iw] / radref[iw];
            }

            break;

        default:

            /* Multi-Scattering with Fixed Model */
            /* --------------------------------- */

            *modmin = aer_opt-1;
            *modmax = aer_opt-1;
            *modrat = 0.0;

            if (*modmin < 0 || *modmin > input->naermodels-1) {
                HY1C_out("Invalid aerosol option: %d\n",*modmin);
                exit(1);
            }

            /* convert input NIR aerosol radiance to relectance */
            rhoa[iwnir_l] = La1[iwnir_l]*radref[iwnir_l];

            /* get aerosol reflectance in visible for fixed model, and epsilon */
            status = fixedaer(sensorID,*modmin,wave,nwave,iwnir_s,iwnir_l,
                    input->aermodels,input->naermodels,
                    solz,senz,phi,wv,rhoa,eps);

            /* convert aerosol relectance to radiance */
            if (status == 0)
                for (iw=0; iw<nwave; iw++)
                    La2[iw]  = rhoa[iw] / radref[iw];
            break;
    }

    /* compute diffuse transmittance through aerosol and Rayleigh */
    if (status == 0 && aer_opt != AERRHNIR && aer_opt != AERRHFRNIR && aer_opt != AERRHSWIR && aer_opt != AERRH) {
        diff_tran(sensorID,wave,nwave,iwnir_l,solz,senz,phi,wv,pr,taur,
                *modmin,*modmax,*modrat,rhoa,taua,t_sol,t_sen);
    }

    /* transfer outputs per wavelength to outputs per band */
    for (iw=0; iw<nwave; iw++) {
        La2_out  [iw] = La2  [iw];
        taua_out [iw] = taua [iw];
        t_sol_out[iw] = t_sol[iw];
        t_sen_out[iw] = t_sen[iw];
    }
    *modmin = *modmin+1;
    *modmax = *modmax+1;

    return (status);
}



/* --------------------------------------------------------------- */
/* get_angstrom.c - compute angstrom coefficient                   */
/*                                                                 */
/* Inputs:                                                         */
/*     l2rec - level-2 structure containing one complete scan      */
/*             after atmospheric correction.                       */
/*     band  = band number 0-7                                     */
/* Outputs:                                                        */
/*     angst - angstrom coefficient                                */
/*                                                                 */
/* Written By: B. A. Franz, SAIC, August 2004                      */
/*                                                                 */
/* --------------------------------------------------------------- */
void get_angstrom(l2str *l2rec, int band, float angst[])
{
    static int   firstCall = 1;
    static int32_t  ib2;
    static float wave2;

    int32_t  ip;
    int32_t  ib1;
    float wave1;
    float aot1;
    float aot2;

    if (firstCall) {
        ib2 = windex(865.0,l2rec->fwave,NBANDS);
        wave2 = l2rec->fwave[ib2];
        firstCall = 0;
    }

    if (band < 0)
        ib1 = windex(443.0,l2rec->fwave,NBANDS);
    else
        ib1 = band;
    wave1 = l2rec->fwave[ib1];

    for (ip=0; ip<l2rec->npix; ip++) {
        aot1 = l2rec->taua[ip*NBANDS+ib1];
        aot2 = l2rec->taua[ip*NBANDS+ib2];
        if (aot1 > 0.0 && aot2 > 0.0)
            angst[ip] = -log(aot1/aot2)/log(wave1/wave2);
        else if (aot1 == 0.0 && aot2 == 0.0)
            angst[ip] = 0.0;
        else {
            angst[ip] = BAD_FLT;
            l2rec->flags[ip] |= PRODFAIL;
        }
    }

    return;
}


/* --------------------------------------------------------------- */
/* get_ms_epsilon.c -                                              */
/* --------------------------------------------------------------- */
void get_ms_epsilon(l2str *l2rec, float eps[])
{
    int32_t  ip, ipb1, ipb2;

    for (ip=0; ip<l2rec->npix; ip++) {
        ipb1 = ip*NBANDS+iwnir_s;
        ipb2 = ip*NBANDS+iwnir_l;
        if (l2rec->La[ipb2] > 0.0) {
            eps[ip] = l2rec->La[ipb1]/l2rec->La[ipb2]
                * l2rec->Fo[iwnir_l]/l2rec->Fo[iwnir_s];
        } else {
            /* "should" indicate ATMFAIL, but just to be safe */
            eps[ip] = BAD_FLT;
            l2rec->flags[ip] |= PRODFAIL;
        }
    }

    return;
}
